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ABSTRACT 


Curvature  affects  the  flow  of  water  in  open  channels  in 
several  ways.  These  effects  include  superelevation, 
secondary  (induced)  flows,  longitudinal  velocity 
redistribution,  and  bed  scour  in  the  case  of  a  mobile 
boundary  such  as  a  natural  stream  or  river.  It  is  important 
for  river  engineers  to  understand  and  be  able  to  estimate 
these  effects. 

This  study  first  presents  and  discusses  the  present 
state  of  knowledge  and  understanding  of  the  problem  of 
curved  channel  flow.  The  flow  mechanics  are  discussed  first 
with  an  emphasis  on  the  prediction  of  specific  flow 
phenomena.  The  state  of  the  art  is  found  to  be  adequate 
with  regard  to  some  effects  but  not  others.  A  body  of 
experimental  data  obtained  under  controlled  conditions  is  an 
important  component  of  the  state  of  knowledge  of  a  field. 
The  existing  curved  channel  experiments  are  tabulated  and 
evaulated  in  the  second  part  of  this  section.  A  discussion 
of  the  relatively  recent  application  of  numerical  modelling 
is  presented.  Aspects  of  the  problem  requiring  further 
research  and  elaboration  are  also  identified. 

A  theoretical  analysis  is  then  performed  based  on  a 
perturbation  of  the  governing  flow  equations.  This  analysis 
is  confined  to  rectangular  channels  but  may  be  extended  to 
an  arbitrary  channel  shape.  The  theoretical  predictions  of 
developed  and  developing  longitudinal  velocity  distributions 
and  maximums  are  found  to  agree  well  with  experimental 


IV 


measurements . 


An  experimental  study  comprising  two  complete  runs  in  a 
single  long  (  270°)  bend  was  also  performed.  Use  of  the 
Laser  Doppler  Anaemometer  to  measure  longitudinal  and 
lateral  velocity  components  allowed  unprecedented  detail  and 
precision.  Analysis  of  the  velocity  data  provided  the 
distributions  of  bed  stress  and  bed  stress  angle  which  are 
important  for  scour  prediction. 

The  purpose  of  this  study  was  to  apply  recent  research 
to  practice  and  to  attempt  to  fill  a  few  of  the  gaps  in  the 
present  state  of  knowledge  of  curved  channel  flows. 
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INTRODUCTION 


Curvature  effects  in  open  channels  affect  many  problems 
in  hydraulic  engineering  such  as  flow  distribution,  bank 
stability,  bed  topography,  resistance,  and  superelevation. 
In  order  to  provide  more  reliable  and  economical  engineering 
works,  these  effects  must  be  taken  into  account. 
Unfortunately,  the  three-dimensional  nature  of  the  problem 
makes  it  very  difficult  to  solve  it  analytically.  Due  to 
the  large  numbers  of  variables  and  geometries  inherently 
possible,  a  purely  empirical  approach  is  not  feasible 
either. 

Recently,  because  of  the  importance  of  the  problem, 
much  research  has  been  conducted  on  flow  in  curved 
channels.  Typically  the  researchers  have  attempted  to  apply 
an  approximate  analysis  and  performed  a  limited  experimental 
study.  The  mechanics  of  the  flow  have  thus  been  revealed 
piecemeal  as  each  study  pointed  to  new  aspects  of  the 
problem.  The  present  knowledge  is  substantial,  but 
fragmented.  The  depth  of  understanding  is  also  uneven  with 
developed  secondary  flow  well  documented  while  main  velocity 
redistribution  remains  difficult  to  predict.  What  knowledge 
that  is  available  is  not  of  a  form  useful  to  the  practicing 
engineer  due  to  the  incomplete  and  fragmentary  (and 
occasionally  contradictory)  nature  of  the  research. 

In  this  study,  it  is  proposed  to  first  review  the 
present  knowledge  of  curved  channel  flows.  Both  the  extent 
and  applicability  of  the  results  are  assessed.  Weak  areas 
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may  thus  be  identified  for  further  study. 

The  subsequent  parts  of  this  study  constitute  a 
contribution  to  one  such  relatively  weak  area;  that  of  the 
developing  flow  situation.  For  simplicity  a  rectangular 
channel  shape  is  used.  An  approximate  anaylsis  is  developed 
and  two  detailed  experiments  are  presented. 


1. 


PART  ONE 


REVIEW  OF  THE  PRESENT  KNOWLEDGE 


1.1  INTRODUCTION 

Almost  all  rivers  and  streams  are  curved  or  meandering 
to  some  degree.  This  curvature  has  a  considerable  effect  on 
the  flow,  bed  topography,  and  stability  of  the  stream. 
Usually,  in  practical  situations,  designers  avoid  curved 
reaches  or  take  account  of  these  effects  with  simple  but 
overconservative  estimates.  It  is  the  purpose  of  this 
section  to  analyse  the  state-of-the-art  knowledge  of  curved 
channel  flows  and  to  obtain  estimates  of  these  effects  such 
that  may  be  used  in  hydraulic  engineering  practice. 

Generally,  the  state  of  present  knowledge  may  be 
divided  into  the  following  three  areas: 

1.  The  flow  mechanics  which  form  the  basis  of  our 
understanding  of  the  problem. 

2.  A  body  of  experimental  information  which  provides 

validation  of  theories  and  guidelines  for 

practical  situations. 

3.  Predictive  models  which  may  be  used  for  detailed 
design  and  analysis. 

Each  of  these  will  be  considered  separately. 

Related  topics  which  are  very  interesting  but  outside 
the  present  scope  are  meander  origin,  development,  and 
migration  as  well  as  dispersion  and  diffusion 
characteristics  of  curved  channels. 
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1.2.  MECHANICS  OF  CURVED  CHANNEL  FLOW 


1.2.1  Introduction 

In  order  to  make  engineering  judgements  and 
calculations,  a  reasonable  understanding  of  the  physics  of 
the  problem  is  required.  At  the  present  time,  some  aspects 
of  flow  in  curved  open  channels  are  well  understood  while 
others  are  not.  The  present  section  will  attempt  to  explain 
the  features  of  curved  channel  flow  in  an  engineering 

manner;  that  is,  to  provide  a  qualitative  description  which 
is  supported  with  straight  forward  formula  derivations  and 
pertinent  experimental  results. 

1.2. 1.1  Co-ordinate  System 

The  co-ordinate  system  used  will  be  based  on  Figure  1, 
where  x  and  u  denote  the  longitudinal  direction  and 

(time-averaged)  velocity,  y  and  v  the  lateral  direction  and 
velocity,  while  z  and  w  denote  the  vertical  direction  and 

velocity.  The  radius  of  curvature  of  the  channel  centreline 
is  R  and  the  width  and  average  depth  of  the  channel  are 

given  by  B  and  H  respectively. 

The  radius  of  curvature  of  any  point  across  the  channel 
is  given  by  r  which  is  related  to  R  and  y  by  the  relation 

r  =  R  +  y  ( 1 ) 

The  elevation  of  the  water  surface  above  H  (referring  to 
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Figure  2)  is  h  while  Hs  is  the  scour  depth  below  the  mean 
bed  elevation.  The  actual  depth  of  flow  is  represented  by  a 
and  is  given  by: 

d=H+Hs+h  (2) 

1.2. 2.1  General  Assumptions 

All  velocities  are  assumed  to  be  time  averaged  and  the 
flow  itself  is  assumed  to  be  steady.  The  depth  averaged 
longitudinal  velocity  is  U  and  the  cross-sectional  average 
is  given  by  Uq. 

In  much  of  what  follows,  the  flow  will  be  assumed  to  be 
"developed",  that  is,  invariant  in  the  longitudinal 
direction.  This  condition  has  also  been  described  as 
" axisymmetr ic"  .  The  requirements  for  developed  flow  are  as 
f ol lows : 

1.  The  radius  of  curvature  of  the  channel  (R)  remains 
constant  (i.e.,  the  bend  is  circular  in  plan). 

2.  The  channel  is  prismatic. 

3.  The  section  under  consideration  is  well  removed 
from  the  entrance  and  exit  (i.e.,  the  bend  is 
long ) . 

These  conditions  are  quite  restrictive  and  are  rarely 
met  in  a  laboratory,  much  less  in  the  field.  The  developed 
flow  assumption  is  useful,  however,  in  isolating  some 
features  of  the  flow  for  better  understanding. 
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The  procedure  in  the  following  will  be  to  attack  the 
individual  phenomena  separately.  This  results  in  the 
simplest  presentation,  but  it  must  be  remembered  that  there 
are  significant  interactions  between  them  and  a  proper 
approach  would  be  to  find  a  solution  for  the  flow  as  a 
whole.  The  relative  depth  of  presentation  of  the  different 
aspects  is  indicative  of  the  state-of-the-art  understanding. 

1.2.2  Superelevation 

The  most  obvious  consequence  of  flow  around  a  curved 
channel  is  superelevation,  that  is,  the  rise  of  the  water 
surface  around  the  outside,  and  the  fall  of  the  water 
surface  around  the  inside  of  the  bend.  The  magnitude  of 
this  water  surface  slope  is  not  large  but  it  has  important 
indirect  consequences. 

The  cause  of  this  phenomena  is,  of  course,  centrifugal 
force.  For  the  water  to  follow  a  circular  path,  a  force 
must  be  applied  in  a  direction  normal  to  the  flow  and 
directed  toward  the  centre  of  curvature.  This  is  provided 
by  the  excess  hydrostatic  pressure  of  an  increased  depth  of 
water  towards  the  outside  of  the  bend. 

1.2. 2.1  Derivation 

To  calculate  the  water  surface  profile,  consider 
Figure  3,  a  free  body  diagram  of  a  water  column  of  unit 
thickness  flowing  concentrically  around  a  bend  with  depth 
averaged  velocity  U.  The  balance  of  forces  acting  on  this 


9 


Figure  3. 


Forces  Acting  on  a  Water  Column 
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fluid  element  is: 


mU2 

r 


F  =  0 


(3) 


where  PQ  and  are  hydrostatic  pressure  forces  and  F  is 

lateral  bed  shear  force.  The  mass  of  the  element,  m,  is 
given  by  the  expression 


m  =  —  d  Ar  (  4  ) 

g 

where  y  is  weight  per  unit  volume  of  the  fluid  and  g  is  the 
acceleration  due  to  gravity.  Substituting  4  and  the 

hydrostatic  forces  into  Equation  3  we  have: 

2 

U  y  jA  r  y  ,  ,  , Ah.  2  y  ,,  Ah. 2^  ~  /c. 

—  ^  dAr  ■-  [-£  (d+— j)  -  j  (d— j)  J  -  F  =  0  (5) 

2 

Simplifying  and  ignoring  terms  of  order  (Ah)  : 

n2  v 

—  1  dAr  -  yd Ah  -  F  =  0  (6) 

r  g 


or : 

Ah  _  U2  F 

Ar  gr  yd Ar 

or  in  the  limit  as  Ar  -*  0: 

dh  _  Ty0 

cfr  gr  dr 


(8) 


11 


where  x  ^  is  the  bed  shear  stress  in  the  lateral 

yo 

direction.  In  many  cases,  it  is  small  and  may  be 
neglected.  As  a  first  approximation  we  may  take  U  equal  to 
Uq  and  r  equal  to  Rq  where  Rq  is  the  radius  of  curvature  of 
a  long  circular  reach. 


dh 

dr 


(9) 


Equation  9  may  be  integrated  by  noting  from  Equation  1: 


dr  =  dy 


(  10) 


to  give: 


h 


0__ 

gR0 


(li) 


At  the  outside  of  the  bend: 


h 


(12) 


and  the  total  superelevation  (difference  between  outside  and 
inside  elevations)  is: 


B 

R 


0 


S 


2 


(13) 
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1.2. 2. 2  Superelevation  Example 


As  an  example  to 

illustrate 

the  size  of 

this 

superelevation,  take  the 

North 

Saskatchewan 

River 

bend 

around  the  Mayfair  golf 

course 

in 

Edmonton 

(data 

from 

Kellerhals,  et  al.,  1972), 

where : 

B  -  140  m 


(14) 


UQ  =  1.25  m/s 


(15) 


Rq  =  750  m 


(16) 


Thus  : 


S  -  30  mm 


(17) 


This  does  not  seem  large,  considering  that  this  particular 
bend  is  of  only  moderate  curvature,  but  the  cross  channel 
water  surface  slope: 


|  =  0.00021 

D 


(18) 


is  of  the  same  order  of  magnitude  as  the  longitudinal  bed 
slope : 


S 


0 


0 . 00035 


(19) 
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The  above  estimation  of  superelevation  has  been  found 
experimentally  to  be  accurate  to  within  about  ten  percent, 
although  it  appears  to  consistently  underestimate  the 
superelevation.  In  some  cases  of  erodable  bed  laboratory 
models,  superelevations  of  up  to  twice  this  value  have  been 
observed  (Yen  and  Yen,  1971). 


1.2. 2. 3  Effect  of  Varying  Longitudinal  Velocity 

To  improve  this  estimate,  one  approach  taken  is  to 
allow  U  to  vary  across  the  channel.  The  problem,  then,  is 
to  specify  the  variation  of  U.  Ippen  and  Drinker  (1962) 
tried  several  variations: 


a.  the  free  vortex;  U 


U 


R 

0  F 


0 


(20) 


Just  past  the  entrance  of  the  bend,  the  velocity  will 
be  distributed  according  to  equation  20.  The  corresponding 
expression  for  superelevation  was  derived  to  be: 


S  = 


B  U0 
R0  29 


1 


(_!_) 

2rq 


b.  the  forced  vortex;  U  =  U 


0  R 


0 


(21) 


(22) 


The  "forced  vortex"  distribution  approximates  the  flow  near 
the  end  of  the  bend  where  the  maximum  velocity  has  shifted 
to  the  outside.  The  superelevation  is  then: 
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1 


2 


(23) 


0 


Another  possibility  is  to  allow  r  to  remain  variable.  This 
results  (assuming  U  =  Uq  again)  in: 


(24) 


Figure  4  is  a  plot  of  Equations  19  ,21,  23  and  24.  The 
different  formulas  diverge  widely  for  large  values  of  b/Rq  , 
but  at  values  typical  of  natural  streams  (B/Rq  <  0.5) 


the  difference  is  very  small 


The  free  vortex  form 


Equation  21,  is  generally  preferred  since  it  predicts  a 
superelevation  slightly  larger  than  the  simple  expression 
given  by  Equation  19.  The  actual  cross-channel  profile  is 
also  best  fitted  by  the  free  vortex  form  ( Ippen  and  Drinker 
1962;  Yen  and  Yen,  1965;  Yen,  1971). 

1.2. 2. 4  Effect  of  Bed  Friction 

It  remains,  however,  to  consider  the  effect  of  the  bed 
friction  on  the  superelevation.  This  bed  friction  acts  as  a 
result  of  the  secondary  circulation  which  is  to  be 
considered  later.  For  now,  we  will  borrow  two  significant 
results  from  this  later  section  that  will  allow  us  to 
estimate  a  correction  to  the  superelevation. 

First,  the  secondary  flow  near  the  bed  is  always 
towards  the  inside  of  the  bend.  This  implies  that  the  bed 


shear  force  acts  in  a  negative  sense  according  to  the  arrow 


Equation 
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Figure  4.  Superelevation  Versus  Curvation  Ratio 
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on  Figure  3.  Thus,  from  Equation  8,  the  lateral  water 
surface  slope  is  always  effectively  increased  by  the  bed 
friction. 

To  estimate  the  size  of  the  correction  let: 


t  A  =  t  Atan<5 
yO  xO 


(  25) 


where  is  longitudinal  bed  stress  and  5  is  the  angle  of 
the  total  bed  stress  vector  from  the  longitudinal 
direction.  Also  define: 


Tx0  =  -(Y/g)  u* 


(26) 


where  u*  is  the  friction  velocity  and 


C2  =  U2/u2 


(27) 


where  Cf  is  a  nondimens ional  Chezy  coefficient 


Cf  =  C//g  ( 28) 

Substitute  the  above  into  Equation  8: 

dh  _  tan  5  (2g) 

Sr  -  qr  g  Q26 

Now  make  use  of  the  second  result  to  be  obtained  later. 


namely: 
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tan6  -  11  — 
r 


(30) 


Equation  29  becomes: 


dh 

dr 


(31) 


Since  typical  values  of  Cf  range  from  about  10  (rough)  to  20 
(smooth),  the  correction  is  about  2.5%  to  10%  and  accounts 
for  some  of  the  observed  error.  This  correction  applied  to 
any  of  Equations  19  to  24  will  result  only  in  the 
magnification  of  the  curves  by  the  factor  (1  +  11/C^).  In 
practice  it  is  difficult  to  apply  this  correction  because 
the  constant,  taken  above  as  11,  actually  is  a  function  of 
Cf  itself  and  has  been  quoted  as  varying  from  about  7  to  13 
( Koch ,  1980 ) . 


1.2. 2. 5  Calculation  of  Superelevation 

In  practical  river  situations  where  it  is  of  interest 
the  superelevation  may  be  calculated  by  the  following 
method : 

1.  Obtain  the  average  stage  elevation  for  the  desired 
discharge  by  standard  methods. 

2.  Obtain  a  cross-sectional  bed  profile  of  the 
channel  (which  should  be  measured  but  may  be 
estimated  using  the  results  of  section  2.6). 
Variation  of  the  resistance  coefficient  across  the 
channel  (e.g.  flood  plains)  should  be  noted  as 
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well. 

3.  Calculate  U  at  a  number  of  evenly  spaced  locations 
across  the  channel.  At  present,  the  best  estimate 
is  given  by  a  resistance  equation,  such  as: 


U  =  Cf  /gdSQ  (32) 

Note  that  the  local  resistance  coefficient  may 
require  a  minor  adjustment  to  give  the  same 
discharge  as  specified. 

4.  Integrate  Equation  24  numerically  across  the 

channel  starting  at  the  inside  bank.  Use  the 
values  of  U,  r  (measured  from  the  centre  of 

curvature)  and  estimated  previously  at  each 

point.  This  gives  the  water  surface  profile 
relative  to  the  elevation  at  the  inside  bank. 

5.  The  constant  of  integration  is  determined  by 

calculating  the  average  elevation  increase  of  the 
profile  from  step  4.  This  value  is  subtracted 

from  the  water  surface  elevations  to  give  h(y). 

6.  It  may  be  necessary  to  repeat  the  process 

(starting  at  step  3)  once  using  the  new  water 
surface  elevations. 


1.2. 2. 6  Developing  Flow  Considerations 

For  the  case  of  developing  flow,  little  is  known.  At  a 
sudden  change  in  curvature,  such  as  at  the  beginning  or  exit 
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of  a  circular  bend,  the  water  surface  slope  has  been  noticed 
to  adjust  over  a  length  about  equal  to  the  width  of  the 
stream: 


lb  *  B 


(33) 


and  roughly  symmetrical  on  each  side  of  the  point  of 
curvature  change.  Over  a  long  bend,  the  magnitude  of 
superelevation  decreases  slightly  as  the  maximum  flow 
concentration  shifts  from  the  inside  to  the  outside  of  the 
bend.  This  corresponds  to  a  gradual  shift  of  superelevation 
from  a  value  given  by  Equation  21  to  that  given  by  Equation 
23. 

The  lateral  inertia  of  the  flow  may  also  have  a 
significant  effect  on  the  transverse  water  surface 
profile.  Yen  and  Yen  (1971)  and  Kalkwijk  and  de  Vriend 
(1980)  both  noticed  a  large  increase  in  superelevation  (up 
to  twice  the  normal  value)  in  experiments  performed  in 
channels  with  a  rigid  but  "natural"  bed.  That  is,  there  was 
a  scour  hole  effect.  Both  experimenters  attributed  this 
effect  mainly  to  lateral  inertia.  The  mechanism  is 
essentially  as  follows.  As  the  flow  proceeds  into  the  bend, 
it  encounters  the  scour  hole  on  the  outside  of  the  bend. 
Here  the  flow  is  faster  and,  therefore,  a  large  portion  of 
the  discharge  shifts  to  the  outside.  This  shifting  envolves 
considerable  inertia  and,  in  order  to  stop  the  lateral  flow, 
an  extra  transverse  water  surface  slope  is  set  up. 
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This  effect  has  been  taken  into  account  only  in  very 
complex  numerical  models  (Leschziner  and  Rodi,  1979;  de 
Vriend,  1981).  As  it  is  expected  that  in  some  cases  the 

effect  is  significant,  especially  in  natural  channels,  more 
research  and  application  is  required. 

1.2.3  Secondary  Flow 

Secondary  circulation,  or  helicoidal  flow  as  it  is 

sometimes  called,  arises  from  the  fact  that  some  of  the 
fluid,  specifically  near  the  bed,  is  moving  slower  than  the 
average  velocity.  Since  it  is  moving  slower,  it  lacks  the 
apparent  centrifugal  force  to  resist  the  imposed  pressure 
gradient  of  the  water  surface  slope.  It  is,  therefore, 

accelerated  inward  until  friction  provides  a  balancing 
force.  Fluid  moving  faster  than  average  (near  the  surface) 
would  experience  the  opposite  effect  and  be  accelerated 

outward.  The  net  result  is  a  circulation  over  the  flow 
cross-section,  superimposed  on  the  main  longitudinal  flow 
(Figure  5).  The  total  streamlines  follow  a  helicoidal  path. 

The  secondary  flow  is  important  because  it  is 
considered  to  be  the  primary  cause  of  river  bend  scour.  It 
also  has  significant  effects  on  the  lateral  diffusion 
characteristics  of  a  stream. 

1.2. 3.1  Analysis  of  Developed  Secondary  Flow  in  a  Wide  Open 
Channel 

The  lateral  force  balance  equation  for  an  element  of 


fluid  is 
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Figure  5. 


Secondary  Circulation 
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u2  _  3 (gh )  9 Tyz  g 

r  3r  3  x  y 


(34) 


where  iyZ  is  the  lateral  shear  stress  acting  on  the  top  and 
bottom  faces  of  the  element.  Equation  34  has  been  derived 
and  used  by  Rosovskii  (1961),  Yen  (1965)  and  Engelund 
(1974).  It  is  directly  analagous  to  Equation  8  except  that 
it  has  been  derived  for  a  fluid  element  of  depth  dz  instead 
of  for  a  column  of  depth  equal  to  the  depth  of  flow  d. 
Alternatively,  Equation  34  may  be  derived  from  the  Reynold’s 
Equation  for  the  radial  direction  by  neglecting  the  inertial 
terms  and  the  small  stress  terms. 

To  proceed  further,  substitute  Equation  8  into  Equation 
34  (noting  that  h  is  not  a  function  of  z  so  the  partial 
derivative  may  be  replaced  by  an  ordinary  derivative): 


u2  _  _  +  g  >  +  3xyz  g 

r  r  y  d  3z  y 


(35) 


A  small  change  has  been  introduced  into  Equation  8.  The 
2 

term  u  ,  defined  by: 


/  u2dy 
0 


(36) 


(which  is  the  average  over  the  depth  of  the  longitudinal 
velocity  squared)  has  replaced  the  depth  averaged  velocity 
squared,  U2 ,  which  is  given  by: 

□2  =  (i  /  u  dy)2 
a  0 


(37) 


f.n»  <£»«£>  n.»  ,u»«)  ,,ne 
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The  difference  between  the  two  is  small,  but  the  former  is 
more  correct.  Rearranging  Equation  35: 


9  T 

g  _ y z 

y  9z 


1 

r 


g  Ty° 

Y  d 


(38) 


Let 


u  =  U  +  u' 


(39) 


where  U  is  the  depth  average  velocity  and  u*  is  the 
deviation  from  the  average.  Substitute  into  Equation  38: 


3  T  _ _ 

^  =  -  (u2  +  2  Uu'  +u'2  -  u2-  2Uu '  -u 1  2 ) 

Y  3z  r 

Simplify  Equation  40  to  obtain: 


g  lyo. 

Y  d 


(40) 


9  T 

g  _ yz 

Y  9z 


-  ( -2Uu ' 
r 


u'  ) 


g  Ty° 

Y  d 


(41) 


To  proceed  further  let: 


g  9v 

-2-  T  =  V  - 

Y  y0  t  9z 


(42) 


where  vt  is  the  eddy  viscosity  and  may  be  a  function  of  z. 
Also  let: 


—  x  n  =  -u2tan5 
Y  y  0 


(43) 


\ 
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where  6  is  as  defined  previously.  Assume: 


u ' 


u*f  x 


v 


t 


U.d  f2(|) 


n  =  z/d 


(44) 


(45) 


(46) 


Substitute  Equations  42  to  46  into  Equation  41: 


3  n 


3v  _  Ud  ,  1  2  .2  r  .  . 

2  3n  r  ^  2fl  cr  (fl  fl  ~  d  tan<5)) 


(47) 


Integrate  twice  with  respect  to  n  to  obtain: 


v  = 


Ud 


/  —  [/{-2f1-^-(f2-f|  -  |tan6) }dn+C1]dn+C2  (48) 


The  constants  of  integration  are  determined  by: 


3v 

9n 


=  0  at  n  =  1 


(49) 


(shear  stress  equals  zero  at  water  surface)  and: 


/  vd  n  =  0 
0 


(50) 


(net  lateral  discharge  equals  zero).  It  remains  to  evaluate 
tan  6.  This  term  is  clearly  of  small  magnitude  as  it 
amounts  to  a  small  correction  to  the  lateral  water  surface 
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slope . 

Yen  (1965) 

and 

Rosovskii  (  1961 ) 

ignore  it 

for 

certain 

cases.  It 

may 

be  evaluated  by 

application 

of 

another  boundary  condition  at  the  bed.  The  most  obvious  is 
s imply : 


De  Vriend  (1976)  modified  this  condition  slightly  to  use 
with  logarithimic  velocity  distributions  and  stated: 


v 


(52) 


where  rig  is  the  height  where  the  longitudinal  velocity  given 
by  the  logarithmic  formula  vanishes. 

Another  possible  bed  condition  is  that  the  angle  of  the 
shear  stress  vector  at  the  bed  should  coincide  with  the 
angle  of  the  velocity  vector  near  the  bed: 


v 

u 


=  tan« 

Tx0 


vt  av/Vy 


(53) 


A  difficulty  with  this  condition  is  that  the  level  of  the 
velocity  vector  must  be  specified.  Rosovskii  (1961)  and 
Engelund  (1974)  used  this  condition.  De  Vriend  (1976) 
reviews  these  conditions. 

In  general  the  results  take  the  form: 


K 


d 

r 


tan  6 


(54) 
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where  K  is  a  function  of  Cf  and  of  the  particular  velocity 
and  eddy  viscosity  distributions.  Values  of  K  ranging  from 
7  to  13  have  been  reported,  with  K  =  11  a  reasonable 
approximation  for  typical  stream  conditions  (Koch,  1980; 
Rosovskii,  1961).  Equation  54  may  be  substituted  into 
Equation  48  which  may  be  written  as: 


v 


Ud 

r 


[I1  +  (T] 
f 


(55) 


where  : 


Ix  (  n) 


( /  -  2f . d  +  a. )  dn  +  a0 
1  ri  1  2 


(56) 


and : 


I2  (n)  =  /  ^  C / {“f i  +  +  K  }  dn  +  b1)  dn  +  b2  (57) 

The  problem  is  thus  solved  when  the  distributions  f^ 
and  f2  are  supplied.  Figure  6  is  a  graph  of  1-^  and  I2  when 
the  logarithimic  distribution  of  u  is  chosen  for  f-^  and  a 
parabolic  distribution  of  eddy  viscosity  is  chosen  for  f2. 
K  is  assumed  to  equal  11.  Von  Karman's  constant  is  taken  as 
0.4. 

In  general,  increasing  roughness  (decreasing  )  leads 
to  smaller  secondary  velocities,  especially  near  the  bed. 
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Figure  6. 


Integral  Functions 
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1.2. 3. 2  Bank  Effects 

At  the  stream  banks,  the  lateral  flow  is  turned 
vertically  to  complete  the  secondary  circulation.  At  the 
outside  bank,  the  flow  is  downward  while  at  the  inside  bank 
it  is  upward.  Calculations  have  been  made  for  this  region 
(Rosovskii,  1961;  de  Vriend,  1973)  and  some  useful 
conclusions  have  been  obtained  (for  rectangular  channels  at 
least )  : 

1.  Wall  effects  (that  is,  appreciable  vertical 
velocity)  extend  about  twice  the  flow  depth  into 
the  channel  from  the  banks. 

2.  The  maximum  vertical  velocities  are  about  the  same 
as  the  maximum  lateral  velocities  in  the  centre  of 
the  channel. 

For  other  than  rectangular  channels,  information  is 
scarce . 


1.2.3. 3  Developing  Secondary  Flow 

As  the  fluid  possesses  inertia,  there  is  a  certain 
length  required  to  set  the  secondary  circulation  into  motion 
after  the  bend  entrance.  If  we  assume  that  the  water 
surface  changes  abruptly  from  flat  to  superelevated,  then 
the  development  length,  or  distance  to  reach  90%  of  its 
developed  value,  has  been  estimated  to  be  (Rosovskii,  1961; 
Nouh  and  Townsend,  1979): 
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L  =  (1.8  to  2.3)  H  C, 
s  t 


(58) 


The  corresponding  development  length  for  superelevation  was 
noted  in  Equation  33.  Thus  the  ratio  is: 


For  channels  with  aspect  ratios  (B/H)  of  20  to  40,  the 
two  lengths  are  of  similar  size.  For  larger  aspect  ratios 
(and  rougher  channels)  the  secondary  flow  development  is 
faster  and  the  overall  development  is  governed  by  the 
superelevation  development.  For  smaller  aspect  ratios  or 
smoother  channels  (typical  of  laboratory  situations)  the 
secondary  flow  development  governs. 

At  the  bend  exit  the  same  situation  applies  in  reverse 
and  the  decay  lengths  are  the  same  as  above. 


1.2.4  Longitudinal  Velocity  Redistribution 

The  longitudinal  component  of  velocity  is  significantly 
affected  by  the  curvature  of  the  stream.  As  the  flow  enters 
the  bend,  the  velocity  on  the  inside  increases  while  the 
velocity  on  the  outside  decreases.  As  the  flow  proceeds 
around  the  bend,  the  lateral  profile  becomes  uniform  again 
and  then  begins  skewing  to  the  outside.  If  the  bend  is  long 
enough,  equilibrium  may  be  reached.  Near  the  exit  of  the 
bend,  a  major  redistribution  toward  the  outer  wall  takes 
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place.  This  exit  redistribution  is  equal  but  opposite  to 
the  entrance  redistribution.  Thus  the  largest  "skewness"  of 
the  lateral  velocity  profile  is  actually  observed  just  past 
the  end  of  the  bend. 

There  are  several  different  mechanisms  causing  this 
phenomenon  and  each  will  be  discussed  separately.  In  the 
following,  the  longitudinal  velocity  referred  to  is  taken  to 
be  the  depth  averaged  velocity,  U,  defined  previously. 


1. 2.4.1  Inertial  Effects 

Over  a  short  distance  near  a  change  in  curvature,  a 
velocity  redistribution  takes  place  that  is  essentially 
independant  of  friction.  This  results  in  the  free  vortex 
velocity  distribution  and  may  be  considered  simply  a 
conservation  of  vorticity  from  the  first  curvature  to  the 
second . 

To  calculate  this  effect,  consider  the  specific  energy 
along  a  streamline: 


E 


(60) 


Point  1  is  in  the  straight  reach  upstream  of  the  bend  and 
point  2  is  in  the  bend.  Differentiating  Equation  60  with 
respect  to  r: 


(61) 


But  we  know  already  from  Equation  8  that: 


dr  rg 


Thus : 


dr  r 


Integrating  (noting  that  U2  =  Uq  at  r  =  Rq ) 


U2  = 


U0R0 


(62) 


(63) 


(64) 


as  expected. 

At  a  bend  exit,  we  may  apply  the  same  analysis.  Point 
3  is  before  the  end  and  point  4  is  after  the  end: 


(65) 


For  simplicity  U3  is  assumed  to  be  uniformly  distributed. 
Thus : 


dh3 

dr 


d 

dr 


4 

2g 


(66) 


Again : 


dh 


3 


gr  gr 


dr 


(67) 
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then : 


dU4  U4 


dr  r 


68) 


and  finally: 


69) 


From  the  point  of  view  of  vorticity,  the  phenomenon  may 


be  explained  as  follows.  Initially  the  flow  is  uniform 


across  the  channel  (neglecting  the  sidewall  regions).  Its 
vertical  vorticity  is  therefore  zero.  After  entering  the 
bend,  the  flow  is  constrained  to  move  basically 
concentrically.  This  may  be  achieved  without  altering  the 
vorticity  only  if  the  lateral  profile  assumes  the  form  of  an 
irrotational  (free)  vortex. 

Friction,  of  course,  acts  to  change  the  vorticity  but 
it  requires  a  significant  length  to  cause  an  appreciable 
change.  In  this  case  it  would  act  to  make  the  profile  more 
uniform.  Since  the  flow  is  moving  concentrically,  a  uniform 
distribution  implies  that  the  flow  contains  vorticity.  This 
is  evident  after  the  end  of  the  bend,  where  the  velocity 
profile  skews  outward. 

As  is  obvious  from  the  above,  the  inertial  velocity 
redistribution  is  closely  associated  with  the 


superelevation. 


We  expect,  and  in  fact  notice,  that  the 
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development  lengths  for  velocity  redistribution  and 
superelevation  are  about  the  same.  This  applies  to  both  the 
entrance  and  exit  of  the  bend. 


1. 2.4.2  Slope  Differential  Effect 

Referring  to  Figure  7,  we  observe  that  for  a  given  arc 
0,  flow  on  the  outside  must  travel  further  than  on  the 
inside.  Presumably,  however,  the  change  in  the  bed 

elevation  will  be  the  same  from  A  to  B.  Or: 


S 


0 


Azb 

AX 


Azb 
r  0 


(70) 


For  the  channel  centreline: 


S 


0c 


Az 


V 


Therefore : 


(71) 


S 


(72) 


Thus  the  slope  is  steeper  on  the  inside  of  the 
channel.  To  estimate  the  effect  on  the  lateral  distribution 
of  longitudinal  velocity,  consider  the  uniform  flow 
approximation : 
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Figure  7. 


Slope  Differential 
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U  =  Cfu  =  Cf /gds0 


Substituting  Equation  72: 


U 


/9dS0cR 


0 


✓F 


(73) 


(74) 


But : 


uc  ■  cf 


(75) 


Thus  : 


(76) 


assuming  d  is  constant. 

Equation  76  indicates  that,  if  the  slope  differential 
were  the  only  important  factor,  then,  after  some  distance, 
the  flow  would  still  be  faster  on  the  inside  but  the 
variation  across  the  channel  would  be  reduced.  This  implies 
that  the  flow  gains  vorticity. 


1.2. 4. 2  Depth  of  Flow  (Vertical  Friction) 

In  a  similar  manner  to  the  above,  we  may  include  the 
effect  of  a  non-uniform  bed  profile.  In  Equation  74,  let  d 
be  a  function  of  r  to  obtain: 


■ 
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(77) 


Thus  the  flow  is  faster  where  the  water  is  deeper.  In 
natural  bends,  "it  is  almost  invariably  observed  that  the 
outside  of  the  bend  is  scoured  deeper  than  the  normal  depth 
while  deposition  occurs  on  the  inside  of  the  bend. 
Therefore,  in  natural  channels  a  strong  skewing  of  velocity 
toward  the  outside  wall  may  take  place. 


1. 2.4.4  Lateral  Friction 

The  ultimate  effect  of  lateral  friction  would  be  to 
cause  the  flow  to  move  as  a  rigid  body  around  the  bend. 
This  would  require  that  the  lateral  profile  be  given  by  the 
forced  vortex  equation: 


(78) 


This  effect,  however,  is  not  very  strong  as  may  be  shown  by 
consideration  of  the  magnitude  of  the  lateral  stress 
gradient  compared  to  the  vertical  stress  gradient.  The 
vertical  stress  gradient  may  be  estimated  by: 


v 


t 


(79) 


while  the  lateral  gradient  is: 


. 
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uo 

v  — 7  ~  n  (80) 

9y  c  ET 

Thus  the  lateral  stress  gradient  over  the  vertical  stress 
gradient  is  proportional  to  (d/b)2  and,  since  aspect  ratios 
of  natural  channels  are  of  the  order  10  to  far  greater,  the 
lateral  friction  cannot  have  much  effect. 

Near  the  banks,  lateral  friction  becomes  much  more 
important,  especially  in  slowing  down  the  initial  fast  flow 
on  the  inside  bank.  Otherwise  it  has  little  but  local 
effect,  although  it  is  useful  to  note  that  it  does  act  in 
such  a  way  as  to  shift  the  location  of  maximum  velocity  away 
from  the  inside  of  the  bend. 

1.2. 4. 5  Lateral  Momentum  Transport  by  Secondary  Flow 

The  secondary  circulation  effect  on  the  longitudinal 
flow  is  the  least  understood,  but  in  some  cases  the  most 
important  mechanism  of  flow  redistribution.  Essentially  the 
secondary  flow  convects  the  faster  moving  upper  layer  flow 
toward  the  outside  while  moving  the  slower  bed  layer  flows 
toward  the  inside.  Of  itself,  it  accomplishes  nothing  since 
the  fluid  is  merely  replacing  other  fluid  of  similar 
momentum.  If  there  is  a  velocity  gradient  across  the 
channel,  then  the  lower  transport  will  act  to  partially 
counteract  the  upper  transport. 

When  the  fluid  nears  the  walls,  however,  the  vertical 
transport  becomes  important  and  the  effects  are 
significant.  At  the  outer  wall,  fast  fluid  (continuously 
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replenished  from  the  centre  of  the  channel)  is  convected 
downward  causing  an  increase  in  mean  velocity,  while  the 
opposite  occurs  on  the  inside.  The  effect,  first  noticed 
near  the  walls,  would  be  convected  again  by  the  secondary 
flow  towards  the  centre.  Kalkwijk  and  de  Vriend  (1980)  took 
this  effect  partially  into  account  while  only  the  numerical 
analyses  of  varying  complexity  by  Rosovskii  (1961), 
Leschziner  and  Rodi  (  1979)  and  de  Vriend  (1981)  take  full 
account. 

At  present,  the  only  results  are  qualitative  and 
suggest  that  this  effect  is  more  important  for  channels  of 
small  aspect  ratio  (narrow  channels)  and  also  more  important 
for  smooth  channels  simply  because  frictional  effects  are 
smaller  by  comparison.  Therefore,  this  effect  would  be  much 
more  important  in  a  laboratory  flume  than  in  a  natural 
stream . 

1. 2.4.6  Developing  Flow 

As  may  be  observed  from  the  above,  there  remains  much 
that  is  not  quantified  regarding  the  problem  of  longitudinal 
velocity  redistribution.  This  makes  the  situation  of 
developing  flow  very  difficult  to  handle  in  a  practical 
manner.  It  has  been  observed  that  a  large  portion  of  every 
bend,  even  in  the  laboratory,  is  in  a  developing  situation 
with  regard  to  the  longitudinal  velocity  distribution.  The 
superelevation  and  secondary  flow  develop  quickly  and  remain 
fairly  stable  but  the  longitudinal  flow  develops  slowly. 


■ 

■ 
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Further,  it  is  the  rare  natural  bend  that  maintains  a 
constant  curvature.  For  these  reasons,  a  model  of 

developing  flow  is  desired.  Short  of  using  a 

three-dimensional  numerical  model,  a  reasonable 
approximation  is  (Engelund,  1974): 

■  -  y  a*  (k>  -  i  (2U'  -  H'>  (81) 

where : 


U' 


U  -  u0 


(82) 


H' 


(83) 


This  model  takes  account  of  some  inertial  and  bed  friction 
effects  for  channels  with  a  lateral  bed  slope.  It  is, 
therefore,  a  good  approximation  for  wide  natural  streams. 
Momentum  transfer  by  the  secondary  flow  is  not  included. 


1.2.5  Bed  Shear  Stress 

Another  important  area  of  concern  is  the  bed  shear 
stress  distribution.  Ippen  and  Drinker  (1962)  and  various 
others  have  studied  this  problem,  experimentally  at  least. 

As  a  first  approximation,  we  may  take: 

Tx0  2  ,  U2 


P 


(84) 


/ 
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and  assume  Cf  to  be  a  constant  or  a  function  of  the  relative 
roughness  (k  /d). 

From  Equation  84  and  what  we  know  from  the  previous 
section  about  the  behaviour  of  U,  we  can  qualitatively 
predict  the  bed  shear  stress  distribution.  The  shear  stress 
would  be  maximum  near  the  inside  wall  at  the  beginning  of 
the  bend,  then  shift  to  the  outside  wall,  and  finally  reach 
a  peak  near  the  outside  wall  just  past  the  end  of  the 
bend.  This  is  in  fact  what  is  observed  (  Ippen  and  Drinker, 
1962)  . 

The  correlation  is  reasonable,  but  consistent 
discrepancies  do  occur.  For  example,  the  stress  near  the 
inside  wall  is  usually  somewhat  lower  than  predicted  while 
the  stress  near  the  outside  wall  is  somewhat  larger. 
Clearly,  this  is  an  effect  of  the  secondary  flow  on  the 
velocity  distributions  themselves.  Only  in  the  complex 
numerical  models  has  this  been  accounted  for. 

The  lateral  bed  stress  arising  from  the  secondary  flow 
itself  was  evaluated  previously  as: 


x  ~  =  x  p.  K  — 
yO  xO  r 


(85) 


The  total  shear  stress  (magnitude)  is  therefore: 


- 2—  2 

T0  =  /xy0  +  Tx0 


=  T  -  /I  +  (  K-)  2 

xO  r 


(86) 


■ 
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Even  if  we  take  K  to  be  11, 
for  any  natural  stream: 


as  suggested  previously, 


T0  ~  Tx0 


(87) 


1.2.6  Wall  Shear  Stress 

Wall  shear  stresses  are  even  less  well  understood, 
despite  their  importance  in  bank  protection  works.  The  only 
approach  available  seems  to  be: 


w  n2 
-  cc  U 

p  w 


(88) 


where  Uw  is  the  velocity  extrapolated  (by 
(as  if  there  were  no  boundary  layer).  For 
or  non-rectangular  channels  even  this 
doubtful . 


eye)  to  the  wall 
cases  of  natural  * 
approach  becomes 


1.2.7  Friction  Coefficient 

The  average  shear  stress  across  the  bed  in  relation  to 
the  overall  friction  coefficient  is  as  follows: 


i 


0 


P 


(89) 


If  we  average  Equation  84  across  the  channel  we  obtain 
(assuming  Cf  is  constant): 


T 


0 


P 


(u2) 


(90) 


42 


Let  (as  in  Equation  82)  : 


U  = 


V 


(91) 


Substitute  Equation  91  into  90: 


T 


0 


p 


(1  +  U'2) 


Simplify  to  obtain: 


(92) 


C  2 

f  (1  +  u,z) 


indicating  that  the 
decrease  (apparently  i 
the  skewness  of  the 
approximated  by: 


(93) 

overall  friction  coefficient  will 
ncreasing  roughness),  depending  upon 
velocity.  Equation  93  may  be 


C 


2 
f  0 


(94) 


which  may  be  used  to  evaluate  the  overall  friction 
coefficient  if  the  distribution  of  U  is  known.  Many 
studies,  both  analytic  and  experimental  have  been  performed 


to 

determine 

the 

excess 

friction 

loss 

of 

a  curved  reach. 

The 

results 

are 

not  conclusive 

and 

in 

some  cases  are 

contradictory 

A 

better 

understanding 

of 

the  mechanics  of 

the  bed  stress  distributions  is  important  for  this  purpose. 
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1.2.8  Bed  Scour  and  Topography 

Since  we  now  have  at  least  a  reasonable  idea  of  the 
flow  field  in  a  bend,  the  next  step  is  to  attempt  to  predict 
the  stable  bed  topography  or  scour  pattern.  Since,  of 
course,  the  bed  topography  affects  the  flow  field  which  in 
turn  affects  the  bed  topography,  we  expect  that  a  successive 
approximation  technique  will  be  required.  The  method 
outlined  below,  which  follows  Engelund  (1974),  will  be  to 
first  estimate  the  bed  profile  in  a  developed  section,  and 
then  estimate  developing  effects  to  apply  a  correction. 


1.2. 8.1  Developed  Bed  Profile 

Under  developed  conditions,  the  net  rate  of  sediment 
transport  across  the  channel  is  zero.  Thus  the  forces 
acting  on  a  sediment  particle  must  balance.  Figure  8  shows 
the  lateral  forces  acting  on  a  sediment  particle  on  the  bed 
of  a  curved  stream.  The  simplest  force  balance  will  be: 


(W  -  Fl)  sin  a  =  Fdr 


(95) 


where  W  is  particle  weight,  FR  is  lift  force,  a  is  lateral 
bed  angle  and  FDR  is  lateral  drag  force.  FDR  is  assumed  to 
be  proportional  to  bed  shear  stress: 


r 


DR 


tan  5 


(96) 


DL 


where  F 


is  longitudinal  drag  force. 


Fdl  may  be  estimated. 


1 
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Figure  8. 


Particle  Equilibrium 


if  the  particle  is  moving  longitudinally  at  a  constant 
speed,  as: 


45 


F 


DL 


(W  -  F  )  cosatan$ 


where  $  is  the  dynamic  friction  angle  (about 
Substitute  Equation  97  into  Equation  96  to  obtain: 


(W  -  F  )  sina  =  (W  -  F  )  cosa  tan$  tan5 

Lj  Li 


or : 


tana  =  tan$  tanS 


Since : 


tana  = 


dd 

dr 


Equation  99  becomes: 


5^.  =  K  tan  $  - 
dr  r 


Solve  the  differential  equation 


(97) 

25°  )  . 

(98) 

(99) 

(100) 

(101) 


d 


C  r 


(  Ktan<$) 


(102) 
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Substitute  y  for  r  and  solve  for  C  to  obtain: 


d 

K 


K  tan$ 

(i  +|> 


The  maximum  depth  including  scour  is: 


m  . 1  _B  K  tan$ 

H  1  2R 


(  103) 


(104) 


For  moderate  curvature  ratios  (B/Rq  <  0.5)  Equation  104  may 
be  approximated  as: 


,  ,  K  .  R 

1  +  2  tan*  R 


(105) 


Taking  K  as  11  and  $  as  25°,  the  above  becomes: 


1  +  2.5 


B 

R 


(  106) 


This  was  derived  assuming  vertical  sides  and  has  been 
found  to  be  an  overestimate.  An  average  of  field 

observations  (Suga,  1963f  cited  by  Ikeda  et  al.,  1981)  gives 
a  constant  of  approximately  1.5.  Figure  9  shows  a  plot  of 
Equations  104  and  106  with  the  modified  constant  against 
observed  values  for  the  North  Saskatchewan  River  (data  from 
Nwachukwu  and  Neill,  1972).  The  average  depth  and  width 
were  evaluated  at  a  stage  corresponding  to  a  discharge  of 


100,000  cfs. 


. 
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Figure  9. 


Maximum  Scour  Versus  Curvature  Ratio 
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Other  attempts  at  predicting  stable  developed  bend 
scour  have  been  made  by  Kikkawa,  et  al,  (1976),  Koch  (1980 
and  Yen  (1970).  These  involve  more  complete  evaluations  of 
the  force  balances  but  result  in  the  necessity  of  numerical 
integration  and  require  more  detailed  information  about  the 
flow  field. 

Zimmermann  and  Kennedy  (1978)  derived  a  lateral  bed 
slope  based  on  a  cross-sectional  torque  balance.  This 
analysis  showed  that  much  steeper  slopes  than  those 
indicated  by  Equation  101  are  possible  for  certain 
conditions.  Odgaard  (1981)  reviews  these  methods  and 
proposes  a  new  one  based  on  an  incipient  motion  criteria. 

To  obtain  the  bed  topography  for  the  entire  bend, 
including  developing  effects,  the  only  approach  taken  so  far 
is  due  to  Engelund  (  1974).  The  procedure  is  to  set  the 
first  estimate  of  bed  topography  equal  to  a  series  of 
developed  profiles  depending  on  local  curvature.  The 
longitudinal  velocity  field  is  then  calculated  using 
equation  81.  The  mean  lateral  velocities  are  calculated 
from  a  continuity  equation.  Finally,  a  sediment 

conservation  equation  is  solved  numerically  and  the 
correction  is  applied  to  the  initial  estimate.  Engelund 
originally  derived  this  method  for  use  on  a  channel  with  a 
sinusoidally  varying  radius  of  curvature  and  obtained 
reasonable  results  compared  with  the  experiments  of  Hooke 
(1974).  The  method  may  be  generalized  to  other  meander 


patterns . 
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Generally,  the  results  for  sinusoidal  channels  indicate 
the  following: 

1.  The  deepest  scour  is  usually  somewhat  downstream 
of  the  point  of  maximum  curvature. 

2.  The  tip  of  the  point  bar  (deposition)  is  usually 
even  further  downstream. 

Both  the  theoretical  and  experimental  work  supports  these 
conclusions.  A  more  detailed  evaluation  of  this  method  is 
included  in  section  1.4.3. 

1.3.  REVIEW  OF  EXPERIMENTAL  INVESTIGATIONS 

Many  researchers  have  undertaken  studies  on  various 
aspects  of  curved  open  channel  flow.  Due  to  the  complexity 
of  the  problem,  and  its  many  manifestations,  the 
experimental  purposes  and  designs  are  as  diverse  as  the 
number  of  experimenters.  Broadly,  the  experiments  may  be 
classified  into  three  categories: 

1.  Fixed  boundary  channels; 

2.  Mobile  boundary  channels;  and 

3.  Natural  streams. 

Yen  (1965)  gives  a  review  of  experimental  work  prior  to 
his  own.  Included  in  this  report  is  an  update  to  that 
review  containing  the  important  experiments  published  up  to 
1980.  Table  1  summarizes  rigid  boundary  experiments. 
Table  2  is  concerned  with  mobile  bed  results  and  Table  3 
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presents  field  observations  of  natural  streams.  An  overall 
evaluation  of  the  experiments  is  also  included. 


Table  1.  Rigid  boundary  curved  channel  experiments. 
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1.3.1  Fixed  Boundary  Channels 

Most  fixed  boundary  channels  were  constructed  for  one 
of  two  reasons.  First,  as  general  experiments  to  understand 
curved  channel  flows,  and  second  to  evaluate  the  bed  stress 
distribution  for  a  particular  situation.  Thus  some 
experiments  include  detailed  information  on  many  aspects  of 
the  flow  while  others  are  very  limited.  Therefore  the 
diversity  in  the  important  geometrical  and  flow  parameters 
is  not  as  large  as  it  appears.  Additionally,  some  of  the 
experiments  have  unrealistic  values  of  some  parameters.  A 
low  aspect  ratio  ( B/H  <  5)  is  the  most  obvious  example, 
although  a  low  ratio  Rg/H  is  unrealistic  as  well. 

Secondary  velocities  have  not  often  been  measured  due 
to  their  relatively  small  size.  When  measured,  the  lateral 
velocities  have  universally  been  determined  by  measuring  the 
total  horizontal  velocity  and  the  angle  of  deviation  from 
the  longitudinal  direction.  Usually  the  angle  is  measured 
by  lining  up  a  protractor-like  device  with  a  short  thread  in 
the  flow.  The  error  in  these  measurements  is  relatively 
large.  Vertical  velocity  measurements  are  practically 
nonexistent. 

The  only  turbulence  measurements  reported  are  those  of 
Yen  (1965),  but  these  were  taken  in  a  closed  air  duct 
model.  Some  measurements,  especially  of  the  turbulent 
Reynold’s  stresses,  would  be  very  useful. 

An  interesting  result  of  these  experiments  is  that  only 
in  the  configuration  of  Kikkawa  et  al.  (1976)  was  a 
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convincing  demonstration  of  a  fully  developed  situation 
achieved.  The  total  angle  turned  by  their  flume  was  over 
300  degrees. 

i 

1.3.2  Mobile  Boundary  Channels 

Most  mobile  boundary  curved  channel  experiments  were 
set  up  to  investigate  the  bed  topography.  Typically  the  bed 
is  initially  flat  and  then  the  experiment  run  until  an 
equilibrium  bed  topopgraphy  is  reached.  The  bed  elevations, 
and  in  some  cases  various  flow  variables,  are  then  measured. 

Many  of  the  same  criticisms  applied  to  the  rigid 
boundary  channels  may  also  be  applied  to  the  mobile 
channel.  Additionally  it  should  be  noted  that  all  these 
experiments  utilizied  rigid  vertical  walls.  Also,  to 
measure  mean  bed  topography,  it  was  often  necessary  to 
smooth  local  bedforms,  which  may  introduce  further 
distortions . 

1.3.3  Natural  Stream  Observations 

Detailed  studies  of  natural  stream  bends  are  notable  by 
their  relative  scarcity.  The  complexity  and  uncontrolled 
nature  of  the  situation  limit  the  utility  of  the  results. 
Measurements  of  this  kind  are  very  important  in  validation 
of  theoretical  or  experimental  conclusions.  It  is,  however, 
difficult  to  obtain  general  results  from  them. 
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1.4.  MATHEMATICAL  MODELS 

To  calculate  maximum  scour,  superelevation,  and 
secondary  flow,  the  developed  flow  assumption  appears  to  be 
adequate,  provided  that  a  representative  curvature  can  be 
obtained.  For  more  detailed  studies  requiring  bed 
topography  and  flow  distribution  over  a  reach  of  a  natural 
or  artificial  channel,  this  simple  approach  is  inadequate. 
A  more  complex  physical  or  mathematical  model  is  required. 
As  curvature  presents  no  special  difficulties  to  a  properly 
scaled  and  constructed  physical  model,  we  will  focus  our 
attention  on  mathematical  models. 

At  present,  no  complete  mathematical  model  is 
available,  but  some  attempts  have  been  made  at  solving 
various  aspects  of  the  problem.  This  section  will  briefly 
consider  these  approaches  in  order  of  decreasing  complexity. 

1.4.1  Three-Dimensional  Models 

A  fully  three-dimensional  model  based  on  the  full 
Reynold's  equations  and  appropriate  boundary  conditions 
would  be  the  most  complete  possible.  Difficulties  arise, 
however,  in  the  modelling  of  the  turbulent  stresses. 
Rodi  (1980)  gives  a  review  of  the  possibilities,  ranging 
from  simple  constant  eddy  viscosity  assumptions  to  full  two 
parameter  turbulence  transport  models.  Based  on  some 
consideration  of  the  developed  flow  analysis  we  conclude 
that  the  actual  turbulence  model  chosen  is  not  that 
significant.  For  example,  the  secondary  circulation  profile 


obtained  by  assuming  a  mixing  length  model  (Rosovskii,  1961) 
does  not  differ  greatly  from  that  obtained  using  an  average 
eddy  viscosity  model  (Engelund,  1974).  Rodi  (1980)  points 
out  this  observation  as  well  when  considering  a  developing 
strongly  curved  flow. 

The  three-dimensional  model  of  Leschizner  and  Rodi 
(1979)  uses  the  k  -  e  turbulence  transport  model.  This 
model,  although  the  most  complete  yet  attempted,  contains 
some  simplifications.  First,  the  variable  water  surface 
location  is  ignored  in  the  calculation.  That  is,  an 
imaginary  frictionless  lid  is  placed  over  the  flow.  The 
water  elevations  are  calculated  later  from  the  pressure 
distributions.  The  effect  of  this  simplification  is  to 
limit  the  model  to  situations  where  the  Froude  number  is 
small . 

A  second  simplification  was  made  with  regard  to 
longitudinal  derivatives  of  the  turbulent  stresses.  These 
were  neglected,  leaving  only  the  pressure  as  an  effective 
downstream- to-upstream  mechanism.  This  allows  a  storage 
efficient  calculation  procedure,  where  the  pressure  field  is 
first  assumed,  and  then  the  calculation  is  performed 
starting  at  the  beginning  and  proceeding  downstream.  The 
scheme  requires  the  storage  of  only  cross-sectional  velocity 
and  turbulence  values,  reducing  computer  memory  requirements 
considerably. 

This  model  achieved  good  agreement  with  an  experiment 
of  Rosovskii  (1961)  where  the  relative  curvature  was  large 
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(B/Rq  =  1)  and  the  Froude  number  was  small  (F  =  0.1).  The 

channel  shape  was  rectangular  with  rigid  sides  and  bed.  As 
yet,  no  three-dimensional  model  includes  sediment  transport 
features . 

1.4.2  Two-Dimensional  Models 

A  two-dimensional  curved  open  channel  flow  model  could 
be  formed  by  assuming  a  hydrostatic  pressure  distribution, 
and  depth  averaging  the  equations  of  motion.  Assumed 
vertical  profiles  could  be  used  to  calculate  integral 
transfer  coefficients.  Such  a  model  has  been  proposed  by 
Kalkwijk  and  de  Vriend  (1980).  The  flow  is  assumed  to  be 
wide  and  shallow.  The  depth  is  allowed  to  vary  gradually 
across  the  channel,  and  the  channel  need  not  be  prismatic. 
The  effect  of  secondary  flow  convection  on  the  main  flow 
distribution  is  included  in  an  approximate  way. 

The  inertia  of  the  lateral  mean  flow  and  circulation 
are  neglected,  which  restricts  the  model  to  gentle 
transitions  of  curvature  and  bed  topography.  This  model  was 
used  to  give  a  good  prediction  of  the  longitudinal  flow 
field  and  superelevation  in  a  rigid  boundary  flume  of 
natural  bed  topography.  Prediction  of  some  natural  channels 
bends  was  generally  successful  with  some  discrepancies  (de 
Vriend  and  Geldof,  1983). 

Koch  and  Flockstra  (1981)  have  introduced  a  model  which 
also  predicts  bed  topography. 
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1.4.3  One-Dimensional  Models 

If  lateral  distributions  may  be  approximated  by  linear 
functions,  the  problem  may  be  simplified  to  a 
one-dimensional  form.  An  example  is  the  model  proposed  by 
Engelund  (1974)  which  was  partially  discussed  previously  in 
Section  1.2. 4. 6.  The  restrictions  on  this  model  are  at 
least  as  severe  as  those  on  either  of  the  more  complex 
models.  Additionally,  the  restriction  of  gentle  curvature 
(B/Rq  <1)  is  required. 

In  the  original  paper,  Engelund  (1974)  provided  a 
solution  for  the  case  of  sinusoidal  variation  of  channel 
curvature.  This  may  be  easily  generalized  to  an  arbitrary 
alignment  by  solving  Equation  69  numerically.  As  an  example 
of  such  implementation,  the  program  CRVCNL  is  presented  in 
the  appendix  (Section  1.7).  As  of  this  writing,  the 
Engelund  model  is  the  only  one  that  will  calculate  bed 
topography  as  well.  The  main  limitation  of  this  calculation 
is  that  it  applies  to  vertical  rigid  walls  requiring  some 
judgement  to  evaluate  the  scour  depths  if  applied  to  natural 
streams.  The  bed  topography  calculation  has  been  included 
in  CRVCNL  as  well. 

1.4.4  Summary 

Mathematical  modelling  of  curved  open  channels  has  only 
recently  been  attempted  with  some  positive  results.  None  of 
the  models  presented  are  truly  of  a  general  purpose  nature, 
and  indeed  most  are  quite  restrictive. 
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1.5.  CONCLUSIONS 

This  section  has  covered  three  aspects  of  the  present 
knowledge  concerning  the  subject  of  flow  in  curved  open 
channels.  The  following  is  a  brief  evaluation  of  the  state 

i 

of  the  art  in  each  of  these  aspects. 

1.5.1  Mechanics  of  Flow  in  Curved  Channels 

Over  the  past  twenty  years  or  so,  much  effort  has  been 
expended  attempting  to  understand  the  mechanics  of  flow  in 
curved  open  channels.  The  phenomena  of  superelevation  and 
secondary  flow  are  well  understood  for  the  idealized 
situation  of  developed  flow  in  a  wide  channel.  Developing 

flow  and  the  associated  redistribution  of  longitudinal 
velocity  are  not  well  understood,  and  no  simple  predictive 

tools  are  available.  Bed  and  wall  stresses  are  likewise 
poorly  understood.  Equalibrium  bed  scour  is  more 

successfully  predicted,  due  mainly  to  its  dependence  more 

upon  the  secondary  flow  than  on  the  longitudinal  flow  and 
shear  stress.  Even  here,  the  success  is  due  to  the 
calibration  of  the  friction  coefficient  $>,  since  the 
constant  K  and  its  variation  with  C^  are  still  the  subject 
of  considerable  debate. 

1.5.2  Experimental  Studies 

Many  experiments  have  been  performed,  but  most  are 
faulted  by  a  lack  of  scope.  The  nature  of  the  problem  is 

such  that  there  is  an  enormous  range  of  geometrical  and  flow 


variations.  With  a  few  notable  exceptions,  the  present  body 
of  experimental  data  provides  many  geometrical  variations, 
but  with  only  a  few  aspects  of  the  flow  actually  measured. 

As  far  as  the  measurements  themselves  are  concerned, 
there  is  a  need  for  accurate  lateral  velocity,  turbulence, 
and  vertical  velocity  measurements,  all  of  which  are 
practically  nonexistent. 

1.5.3  Mathematical  Models 

Mathematical  modelling  of  curved  channel  flows  is  only 
a  relatively  recent  development.  At  present  there  have  been 
some  promising  approaches  taken.  Much  more  effort  is 
required,  however,  to  extend  these  models  to  practical 
situations  where  arbitrary  geometries,  flow  conditions  and 
sediment  transport  are  important. 


2. 


PART  TWO 


THEORETICAL  ANALYSIS 


2.1  INTRODUCTION 

Developing  flow  in  a  curved  open  channel  is  clearly  a 
complex  three-dimensional  problem.  The  review  of  the 
present  state  of  knowledge  indicated  that  this  problem  is 
however,  important  and  not  well  understood.  It  is  the 
purpose  of  the  present  theoretical  investigation  to  develop 
analyses  which  will  shed  some  light  on  this  difficult  but 
fascinating  problem. 

The  procedure  to  be  followed  will  be  to  first  develop  a 
general-purpose  co-ordinate  system  appropriate  for  a 
meandering  stream.  The  equations  of  motions  are  then 
written  in  this  system  and  non-dimensionalized  using  typical 
scales  for  the  channel  geometry  and  flow.  Assuming 
relatively  mild  curvature,  a  perturbation  analysis  is 
performed  which  results  in  most  of  the  essential  features  of 
the  developing  flow  being  represented  in  the  first  order 
approximation.  An  analysis  is  performed  by  forming 
vorticity-like  equations  to  predict  secondary  flow, 
longitudinal  velocity  redistribution  for  developed  flow  in  a 
bend.  A  simplified  analysis  of  developing  flow  is  also 
done.  Preliminary  comparison  with  experiment  is  also  made. 
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2.2  GOVERNING  EQUATIONS 

2.2.1  Co-ordinate  System 

The  nomenclature  and  co-ordinate  system  to  be  used  in 
the  following  is  the  same  as  was  utilized  in  the  review 
section.  The  co-ordinate  system  will,  however,  be 

generalized  and  formalized  in  a  form  appropriate  for  an 
arbitrarily  specified  meandering  channel.  The  co-ordinate 
directions  are  specified  as  follows. 

1.  The  x  co-ordinate  and  associated  velocity  u  are  in 
a  direction  tangential  to  the  channel  centreline. 
The  x  axis  is,  in  fact,  the  channel  centreline. 
The  origin  may  be  established  arbitrarily.  The 
positive  x-direction  is  downstream. 

2.  The  y  co-ordinate  and  velocity  v  are  perpendicular 
to  the  x  axis  in  plan.  The  y  coordinate  is 
defined  to  be  positive  to  the  left  of  an  observer 
standing  on  the  channel  centreline  facing 
downstream. 

The  z  co-ordinate  and  w  velocity  originate  at  the 
bed  of  the  channel  and  are  positive  upwards. 
Strictly,  the  z  direction  is  perpendicular  to  the 
bed  but  may  be  taken  as  vertical  in  most  cases. 


3. 


For  this  co-ordinate  system  the  metrical  coefficients 
(Nash  and  Patel,  1972)  are: 


_  R ( x)  +  y 

1  RTx 


(108) 


h2  =  1  U09) 

03  =  1  (110) 


R(x)  is  positive  if  the  center  of  curvature  at  that  value  of 
x  is  to  the  right  of  an  observer  facing  downstream. 

This  system  may  be  shown  to  be  triply  orthogonal  by  the 
application  of  the  Lame  conditions  (Nash  and  Patel,  1972). 
As  it  stands,  however,  the  above  system  does  not  establish  a 
unique  relation  between  these  co-ordinates  and  for  example, 
Cartesian  co-ordinates.  In  other  words,  a  given  point  in 
space  may  be  described  by  more  than  one  set  of  co-ordinates 
in  this  system.  In  such  cases,  the  extra  condition  required 
for  uniqueness  will  be  to  choose  that  set  for  which  the 
absolute  value  of  the  y  co-ordinate  is  the  smallest. 

For  use  in  practical  situations  this  system  requires 
the  measurement  of  longitudinal  distances  along  the  channel 
centerline  which  is  the  normal  procedure.  The  variation  of 
radius  of  curvature  must  be  evaluated  by  analytical, 
numerical  or  graphical  means.  An  example  of  a  numerical 
procedure  is  given  as  part  of  the  demonstration  program  in 


the  review  section. 
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2.2.2  Equations  of  Motion 

In  terms  of  the  above  coordinate  system,  the  Reynold's 
Equations  of  motion  for  the  steady  flow  case  with  negligable 
viscous  effects  may  be  written  as: 

Ru  _3u  +  -_8u  -3u  uv  1  _R_  d£ 

R+y  3x  3y  3z  R+y  p  R+y  3x 


R  3u ' 2  aiTV1 

R+y  8x  3y 


au '  w 1 
3z 
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3u  '  v 1 
R+y  ' 


(111) 
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Ru  3w  —  j3w  —  _3w  _1  j3p 

R+y  9z  ay  3z  p3z 


R  3u  1  w  1  +  3v  '  w  1  +  3w  '  2  +  v 1  w  1  \ 
R+y  3x  3y  3z  R+y  ' 


(113) 


where  p  is  the  pieziometric  pressure.  The  overbars  indicate 
time  averaged  quantities  and  the  primes  indicate  fluctuating 
velocity  components. 


bW 


The  corresponding 

continuity 

equation  is: 

R  du  +  8v  3w 

v 

R+y  9x  3y  3z 

+  R+y  "  0 

(114) 

Before  proceeding 

further 

it  is  desirable  to 

non-dimensional ize  Equations  111 

-  114.  Let : 

x*  =  x/b 

(115) 

y*  =  y/b 

(116) 

z*  =  z/dg 

(117) 

u*  =  u/Uq 

(  118) 

v*  =  v/Uq 

(119) 

w  *  =  w/Uq 

(120) 

p*  =  p/ydQ 

(121) 

R*  =  R/Rq 

(  122) 

where  b,  dQ,  Uq  and  Rq 

are  representative  values  (average  or 

maximum  (minimum  in 

the  case 

of  Rg ) )  .  The  following 

non-dimensional  groupings  may  be 

formed  from  these  scales: 

a  =  b/R0 

(123) 

3  =  b/d 0 

(124) 

fR2  =  v02/9d0 

(  125) 

Equations  115-125  may  be  substituted  into  Equations 


111-114  to  give: 
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For  convenience,  from  this  point  on,  all  starred  quantities 
will  be  represented  by  their  unstarred  equivalents. 

The  flow  equations  have  now  been  scaled  such  that  all 
variable  terms  are  of  order  of  magnitude  one  or  less.  A 
preliminary  estimation  of  the  size  of  each  term  is  indicated 
by  the  size  of  it’s  coefficient.  The  three  parameters 
involved  are  a,  3/  and  FR2.  FR2  is  a  typical  Froude  number 


and  may  vary  from  near 

zero 

to 

above  one.  The  parameter  3 

is  usually  referred 

to 

as 

the 

channel  aspect  ratio. 

For 

natural  channels  it 

is 

usually 

vary  large  (~50),  while 

for 

laboratory  channels 

it 

may 

be 

an  order  of  magnitude 

less 

(~10).  Some  experiments  have  been  performed  with  values  of 
3  about  1.  Finally,  a  is  a  measure  of  the  relative 
curvature  of  the  channel.  From  a  maximum  possible  value  of 
one  to  a  minimum  of  zero  it  may  be  the  best  indication  of 
curvature  effects.  Typical  channel  bends  would  have  a  less 
than  0.5.  Average  bends  (Yen,  1965)  are  typically  about 
0.1. 


2.2.3  Perturbation  Expansion 

The  system  of  Equations  125-129  is  not  solvable  by 
present  methods.  In  order  to  derive  useful  results,  it  is 
necessary  to  utilize  approximate  analysis  methods.  A 
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productive  technique  is  that  of  perturbation  analysis  (Lin 
and  Segal,  1972).  This  method  has  been  applied  to  the 
curved  channel  flow  problem  in  the  past  to  predict  secondary 
flow  (Rosovskii,  1961;  Yen,  1965:  de  Vriend,  1976  for 
example).  The  parameter  most  often  chosen  for  the  expansion 
is  a/3  which  certainly  is  very  small.  Since  the  horizontal 
dimensions  were  scaled  with  Rq,  this  resulted  in  drastic 
simplification  of  the  equations  and  postponement  of  any 
consideration  of  developing  flow  to  the  second  order 
approximation  at  least. 

In  the  present  study,  a  different  approach  is 
followed.  The  horizontal  dimensions  are  scaled  with  the 
channel  half-width  which  increases  the  relative  importance 
of  the  differential  terms  involved.  The  perturbation 
expansion  parameter  is  chosen  to  be  a  which  is  larger  than  3 
but  allows  more  information  to  be  gained  from  a  first  order 
approximation.  Still  a  is  small  enough  to  allow  reasonably 
accurate  calculations  for  all  but  the  sharpest  of  bends. 

The  procedure  is  to  expand  all  of  the  unknown  variables 
in  a  power  series  of  a,  which  is  taken  to  be  the  entire 
dependance  of  the  variable  on  a.  For  example: 

—  2  —  2  -  2 
u  (x,y,z, a, 3/Fr  )  =  uQ  (x,y,z,3,FR)  +  au  (x ,y , z , 3 , FR) 


(130) 


Such  series  may  be  formed  for  u,  v,  w,  p,  u 
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v 
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2  _  _ 

w'  ,  u ' w ' ,  and  v ' w ' 

126-129.  In  addition, 


and  substituted  into  Equations 
the  following  approximation  is  used. 
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After  substitution,  straightforward  but  considerable 
algebra  yields  sets  of  equations  for  each  power  of  a.  The 
first  set,  the  zero  order  set,  is: 
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The  first  order  (a1)  set  is: 
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The  zero  order  set  may  be  recognized  as  the  Reynold's 


Equations  in  Cartesian  co-ordinates.  This  indicates  that  to 
zero  order  in  a,  curvature  has  no  effect  on  channel  flow. 
The  zero  order  quantities  may  be  calculated  as  for  straight 
channel  flow.  This  may  be  very  difficult  but  for  the 
present  assume  that  it  can  be  done. 

The  first  order  set  of  equations  look  longer  and  more 
difficult  but  in  fact  are  more  tractable  than  the  zero  order 
set  because  this  set  is  linear  in  the  unknown  first  order 
variables.  The  zero  order  functions  appear  here  as  known 
functions. 

A  particular  simple  situation  of  practical  interest  is 
that  of  uniform  flow  in  a  prismatic  channel.  If  secondary 
flows  due  to  turbulence  are  neglected,  the  solution  to  the 
zero  order  set  is: 
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where 

Sq  is  channel  bed  slope  (143) 

and  all  longitudinal  derivatives  of  the  zero  order 
turbulence  terms  vanish.  For  such  a  case  Equations  136  to 


139  reduce  to: 
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It  may  be  argued  further  that  Equations  144  to  147 


would  also  apply  to  situations  of  gradually  varied 
nonuniform  flow  and  even  to  rapidly  varied  flow  where  the 
changes  are  of  the  order  of  a  or  less.  The  reason  is  that 
these  situations  may  be  constructed  as  a  power  expansion  in 
two  or  more  variables,  for  example: 
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Keeping  only  linear  order  terms  will  give  the  same  two  sets 
of  equations  as  above  plus  an  additional  set  for  each 
additional  variable.  The  different  effects  may  then  be 
superimposed  to  give  the  final  approximate  result. 

There  are,  of  course,  other  ways  of  performing  the 
perturbation  expansions  which  will  result  in  slightly  or 
greatly  different  sets  of  equations.  In  such  a  large  and 
complex  set  of  equations  assumptions  and  approximations  can 
have  unforeseen  consequences. 

Such  a  situation  arises  in  this  particular  set  of  first 
order  equations.  No  assumption  was  made  about  the  size  of 
the  Froude  number  but  the  analysis  is  already  restricted  to 
situations  where  it  is  small.  To  see  this,  consider  the 
situation  in  the  immediate  vicinity  of  a  sharp  change  in 
channel  curvature.  In  such  a  short  location,  friction  may 
be  ignored,  the  channel  may  be  considered  horizontal,  the 
pressure  distribution  hydrostatic  and  zero  order  quantities 
equal  to  a  constant.  Equations  144  to  147  reduce  to: 
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where  u^  and  are  now  depth  averaged  quantities.  These 
three  equations  may  be  solved  for  h^  to  give: 

a2h  a2h 

- 7  +  - i  =  o  (152) 

ax  3y 


Treating  the  problem  a  little  differently  and  depth 
averaging  before  carrying  out  the  perturbation  expansion 
(applying  the  same  set  of  assumptions;  see  Appendix  A) 
results  in: 


< i-fr) 


(153) 


Equation  152  is  ellipitic  for  all  Froude  numbers  while 
Equation  153  becomes  hyperbolic  for  supercritical  flow. 
Equation  153  correctly  predicts  (to  a  linear  approximation) 
the  standing  cross-channel  waves  that  occur  at  a  bend 
entrance  for  supercritical  flow.  For  small  Froude  numbers 
the  two  equations  are  almost  equivalent.  Thus  the  set  144 
to  147  is  limited  to  small  Froude  numbers. 

The  derivation  and  solution  of  Equation  152  is  given  in 
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Appendix  A.  The  results  indicate  that  unless  the  Froude 
number  is  very  close  to  one,  the  effects  are  small.  The 
solution  itself  is  useful  as  it  also  shows  that  appropriate 
length  scale  is  the  channel  half  width  as  was  assumed  in 
this  derivation. 

2.4  TURBULENCE  CLOSURE 

To  proceed  with  calculations  based  on  the  approximate 
set  of  Equations  144  to  147  it  is  necessary  to  relate  the 
turbulence  terms  to  the  mean  flow  quantities.  This  problem 
is  universal  to  turbulent  flow  calculations.  For  the  curved 
channel  flow  situation,  there  are  no  measurements  available 
to  guide  our  approximation.  In  the  review  section,  however, 
it  was  noted  that  the  precise  form  of  the  turbulence  model 
used  was  not  very  important  in  the  calculation  of  the 
'developed'  secondary  flow.  Rodi  (1980)  indicated  the  same 
is  true  for  developing  flow  calculations. 

For  simplicity  the  closure  used  will  simply  be  a 
constant  but  nonisotropic  eddy  viscosity  model  using  a 
straight  wide  rectangular  channel  value.  In  comparsion  to 
other  approximations  made  to  this  point  this  assumption 
should  prove  reasonable.  The  form  for  this  model  is: 
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written  in  tensor  notation  (dimensional  quantities).  A 
proper  non-dimens ional i zat ion  of  is: 


Equation  154  takes  the  following  non-dimensional  form: 


(156) 


This  formulation  assumes  that  the  effective  turbulent 
eddy  viscosity  is  isotropic.  Studies  of  mixing  (Fischer,  et 
al,  1978)  have  shown  however  that  the  lateral  eddy 
diffusivity  (as  distinct  from  eddy  viscosity)  is  at  least 
two  to  three  times  the  vertical  depth  averaged  value.  For 
the  following  a  value  of  three  is  applied  to  the  lateral 
direction  in  the  form  of  a  constant,  vR. 

Equation  156  may  be  substituted  into  Equations  144  to 
147  (stars  dropped)  to  give: 
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where  the  operator: 


V 


2 


+ 


For  the  presant  purpose : 


(160) 


v  =  0.077 


(161) 


From  this  point  on  the  overbars  will  be  dropped. 


2.5  VORTICITY  EQUATIONS 

The  set  of  perturbation  Equations,  Equations  158  to  160 
and  147  can  be  further  reduced  by  eliminating  the  pressure 
by  cross-differentiation.  This  results  in  equations  for  a 
vorticity-like  combination.  The  first  is  derived  by 
differentiating  Equation  157  with  respect  to  y  and  Equation 


158  with  respect  to  x  and  then  taking  the  difference.  After 
some  algebra  the  following  equation  is  obtained: 
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The  term: 


=  (- 


3  v . 
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3y 


(163) 


has  the  appearance  of  vorticity  although  in  this  co-ordinate 
system  it  is  not  equal  to  the  vorticity.  The  form  is 
convenient  however  and  will  be  utilized.  Equation  162 
becomes : 
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The  other  two  "vorticity"  equations  are: 
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where : 
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where : 
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(168) 


These  equations  describe  the  transport  of  perturbation 

"vorticit ies" .  The  vorticities  have  direct  physical 

interpretations  if  the  major  components  are  studied.  The 

longitudinal  vorticity  E,  is  clearly  related  to  the  secondary 

av1 

circulation  through  the  8  — —  term. 
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3u 


3z 


The  lateral  vorticity  through  it's  dependance  on 
,  will  give  information  on  the  distribution  of  bed 


stress.  Finally,  the  vertical  vorticity  n  is  related  to  the 

skewing  of  the  longitudinal  velocity  distriubtion  because  it 

3u, 


contains  a 


3y 


term. 


To  facilitate  discussion  and  highlight  the  important 
features  of  these  vorticity  transport  equations  consider  a 
situation  where  Uq  is  a  function  of  z  only.  This  would 
apply  in  the  central  portion  of  a  wide  channel.  For  this 
situation  Equations  164,  165  and  167  reduce  to: 
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Equation  169  describes  the  transport  of  vertical 
vorticity  which  is  related  to  the  lateral  distribtuion  of 
longitudinal  velocity.  The  rate  of  change  of  this  quantity 
along  the  channel  is  equal  to  the  four  terms  on  the 
right-hand-side  of  Equation  169.  The  first  is  a  turbulent 
diffusion  term  which  would  tend  to  reduce  n  to  0 .  The  third 


term  accounts  for  the  slope  differential  effect  (shorter 
path  around  inside  of  bend) .  The  second  term  provides  the 
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secondary  flow  effect  on  the  longitudinal  velocity.  It  is 
interesting  that  this  effect  is  accounted  for  by  the 
vertical  velocity  distribution.  The  lateral  component 
effect  depends  on  the  lateral  distribution  of  Uq  whereas  the 
vertical  component  effect  depends  on  the  vertical 
distribution  of  Uq  .  The  relative  size  of  this  term  may  be 
estimated  from  previous  studies.  Typically  (Yen,  1965): 


w 
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where  K  =  4 . 


Equation  173  can  be  written  as: 


Therefore : 


(  174) 


(175) 


All  other  terms  are  of  order  one  so  that  for  a  channel  with 


8  -  8f  a  -  0.15  (as  in  the  experiments  to  be  presented 
later) : 
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The  slope  differential  effect  on  the  other  hand  ( FR  *  0.5, 

S0  -  0.001)  : 

6  SQ 

— r —  “  0.03  (177) 

F  R 
R 


And  the  largest  of  the  diffusion  terms  (C^  -  20)  : 


8  v  3u3 

- - ^ -  -  0.05  (178) 

~f  3z  3y 


This  secondary  flow  effect  may  thus  explain  the  skewing  of 
the  velocity  profile  toward  the  outer  bank  even  in 
rectangular  channels. 

The  second  vorticity  transport  equation,  Equation  170, 
concerns  the  longitudinal  component  of  vorticity.  This  is 
associated  with  the  secondary  flow  itself.  By  neglecting 
longitudinal  derivatives,  taking  the  turbulent  viscosity 
ratio  to  be  one,  and  using  the  continuity  equation  to 
provide  a  stream  function,  Equation  170  may  be  written  as: 
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(  179) 


where  : 
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(180) 


This  relation  has  been  derived  and  used  by  Rosovskii  (1961) 
and  others  to  calculate  secondary  flow.  The  forcing  term  is 
related  to  the  radius  of  curvature  and  the  vertical  gradient 
of  longitudinal  velocity.  The  term  including  3v^/3x  in 
Equation  170  can  be  shown  to  represent  the  local  streamline 
curvature,  which  is  especially  important  when  the  boundary 
curvature  changes  suddenly. 

The  last  vorticity  transport  equation,  Equation  171, 
characterizes  the  bed  stress  distribution.  The  last  term 
represents  a  convection  effect.  The  second  last  term  is 
clearly  a  vortex  stretching  effect.  z,  will  be  increased  at 
the  outer  bank  by  this  term  and  decreased  near  the  inner 
bank . 

The  other  terms  involving  Su^/Sy  may  be  shown  to  have 
similar  interpretations.  Overall,  it  appears  that  the 
secondary  flow  currents  interacting  with  the  main  flow 
gradients  have  many  interesting  effects  and  seem  to  account 
for  many  observed  phenomena. 
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2.6  DEVELOPED  FLOW 
2.6.1  Formulation 

Although  the  emphasis  of  this  study  is  on  developing 
flow,  it  is  useful  to  evaluate  the  model  for  an  asymtotic, 
developed  condition.  In  particular,  the  prediction  of 
longitudinal  velocity  distribution  is  of  interest.  The 
relevant  equations  may  be  formed  from  Equations  163  to  168 
by  ignoring  any  terms  involving  a  longitudinal  derivative. 
R  may  be  taken  as  1.0.  These  are: 
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Two  of  these  equations  are  redundant  so  Equations  183 
and  187  will  not  be  considered  further. 

The  calculation  procedure  will  be  to  solve  for  £  from 
Equation  182.  The  secondary  circulation  velocities  may  then 
be  solved  by  using  Equations  184  and  186  as  follows.  Define 
a  stream  function  such  that: 


(188) 


(189) 


This  stream  function  automatically  satisfies  Equation 
184  and  may  be  introduced  into  Equation  186  to  produce: 


(190) 
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Given  E,  from  Equation  182,  Equation  190  may  be  solved  for 
and  then  Equations  188  and  189  give  v1  and  w1#  Equation  181 
then  gives  r\  and  U}  may  then  be  obtained  by  integrating 
Equation  185. 

Boundary  conditions,  especially  at  the  solid 

boundaries,  are  complex  and  controversial.  Some 

simplification  may  be  obtained  by  considering  the  symmetries 
of  the  problem.  If  uq  is  symmetric  about  y  =  0  then  E,  and  \p 
will  also  be  symmetric.  The  vertical  vorticity  n  will  also 
be  symmetric.  The  calculations  may  be  done  then  only  for 
0  <  y  <  1.0.  On  the  virtual  boundary  created  at  y  =  0,  the 
normal  derivative  of  these  quantities  is  zero.  In  addition, 
u^  is  zero  here.  On  the  other  boundaries  \p  may  be  set 

arbitrarily  to  zero,  completing  this  particular  boundary 
value  problem. 

At  the  free  surface,  the  shear  stress  is  approximately 
zero  resulting  in  E,  and  dr\/dz  being  equal  to  zero  at 
z  =  1.  Finally,  boundary  conditions  for  E,  and  n  are 
required  at  z  =  0  and  y  =  1.  As  was  discussed  in  the  review 
section  this  condition  is  somewhat  dependent  on  the 
particular  form  assumed  for  Uq .  Rosovskii  (1961)  and  Ikeda 
(1978)  have  suggested  a  particularly  simple  form  for  E, , 
namely  that  it  will  be  zero  at  the  solid  boundaries.  This 
implies  that  the  shear  stress  is  also  zero  at  these 
boundaries.  It  was  argued  however,  that  this  approximation 
only  affected  the  secondary  flow  distribution  in  the 
immediate  vicinity  of  the  boundary.  This  assumption  makes 
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analytical  solutions,  at  least  for  the  secondary  circulation 
(vj  and  )  possible. 

For  the  present  study,  however,  a  numerical  approach 
was  chosen.  This  allowed  the  most  flexibility  with  regard 
to  distributions  of  Uq  as  well  as  inclusion  of  more 
realistic  boundary  conditions.  Once  a  program  is  written, 
numerical  "experiments"  could  be  performed  to  show  the 
dependence  of  the  perturbation  velocities  on  the  various 
parameters. 

It  is  still  necessary  to  provide  an  initial  velocity 

^u0)  distribution  and  the  remaining  boundary  conditions.  In 
the  interest  of  simplicity,  the  slip  velocity  formulation  of 

Engelund  (1974)  is  a  reasonable  approximation.  Any  other 

prescribed  distribution  could  be  used  but  this  method  is 

consistent  with  the  constant  eddy  viscosity  approximation 

and  leads  to  reasonably  simple  boundary  conditions. 

The  initial  velocity  distribution  may  be  calculated 

from  : 
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with  boundary  conditions: 
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where  is  a  constant  to  allow  the  distribution  to  be  tuned 
to  fit  measured  distributions.  The  constant  Ug^  is  given 
by : 
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This  formulation  has  the  advantage  that  the  same  method 
may  be  used  to  solve  for  the  initial  velocity  distribution 
as  for  the  perturbation  velocities.  The  velocities 
calculated  by  Equations  191-196  must  however  be  normalized 
such  that  the  mean  value  is  1.0. 

The  boundary  conditions  for  the  perturbation  velocities 
may  be  formulated  (Engelund,  1974)  by  the  requirement  that 
the  ratio  of  boundary  perturbation  velocity  to  boundary 
perturbation  stress  be  the  same  as  the  Ug  boundary  velocity 
(slip  velocity)  to  stress  ratio.  In  dimensional  terms  (for 
the  lateral  velocity)  : 
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In  nondimensional  form 
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which  may  be  written  as 
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Similarily : 
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on  y  =  1 
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For  u the  same  condition  gives: 
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which  may  be  converted  into  a  boundary  condition  in  n  by 
differentiation: 
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The  remaining  condition  is: 
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consistent  with  Equation  195  for  the  main  velocity. 

Boundary  condition  202  may  be  used  directly  but  the 
other  three  require  advance  knowledge  of  u^ ,  v-^  and  w^ .  An 
iterative  procedure  must  then  be  used.  As  indicated 
previously  the  effect  of  assuming  zero  boundary  values  does 
not  have  a  large  effect  on  the  overall  solution  so 
convergence  is  very  rapid. 


2.6.2  Solution 

As  indicated  previously  the  above  problem  was  solved 
numerically  for  a  range  of  the  independent  parameters  3,  Cf 
and  n.  A  successive-over-relaxation  iterative  procedure  was 
used  to  solve  the  Poisson-like  equations.  This  procedure  is 
very  simple  to  code,  is  the  most  economical  of  storage 
space,  provides  rapid  convergence  and  works  well  with  the 
iterative  boundary  conditions  imposed.  The  last  point  is 
due  to  the  fact  that  the  latest  approximation  is  used  as  an 
initial  guess  for  the  iteration.  Thus  as  the  boundary  value 
approaches  it's  final  value,  fewer  and  fewer  iterations  are 
required  to  obtain  a  new  value. 

The  FORTRAN  source  code  for  this  problem  is  included  in 
Appendix  B.  Typically,  a  run  took'  about  three  minutes  to 
complete  on  the  University  of  Alberta  Computing  Services 
API  6  4  processor  using  a  fifty  by  fifty  finite  difference 
grid.  Figures  10  to  14  display  the  output  of  the  various 


stages  of  the  calculations.  Figure  15  shows  a  measured 
final  velocity  profile  for  similar  conditions  (from  part  3 
of  this  study ) . 

The  effect  of  vR  and  Cj-,  on  the  output  was  investigated 
and  it  was  found  that  -  0.3  gave  a  reasonable 
approximation  for  Ug  while  vR  =  3.0  gave  reasonable 
results.  Increasing  caused  the  final  perturbation 
longitudinal  velocity  profile  to  be  more  evenly  distributed 
in  the  y  direction  as  is  shown  in  Figure  16.  Decreasing  vR 
results  in  much  smaller  longitudinal  velocity  perturbations 
as  shown  in  Figure  17. 

Using  the  values  quoted  above  the  model  was  run  a 
number  of  times  to  investigate  the  effect  of  the  friction 
coefficient  and  the  channel  aspect  ratio.  The  results  are 
sumarized  in  Figure  18  where  the  depth  averaged  longitudinal 
perturbation  velocities  at  y  =  0.9  are  plotted.  This  graph 
shows  a  strong  dependence  on  aspect  ratio  and  a  weaker 
dependence  on  friction  coefficient. 

Also  investigated  was  the  importance  of  the  various 

forcing  terms  in  Equation  181.  The  overriding  term  was  the 

3w-^ 

one  involving  - — -.  The  rest  of  the  terms  had  the  effect 
of  smoothing  the  lateral  distribution  of  longitudinal 
perturbation  velocity.  Figure  19  shows  the  depth  averaged 
vertical  velocity  at  y  =  1.0  and  while  the  magnitudes  are 
not  large,  the  correspondence  with  Figure  18  is  remarkable. 
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Figure  14.  Total  Longitudinal  Velocity  Contours 

(3  =  8.7,  Cf  =  16.0,  a  =  0 .146) 
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Measured  Longitudinal  Velocity  Contours 
(u/^Of  3  =  8.7,.  C  f  —  16.0.  a  =  0.146) 
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Effect  of  Cb  on  the  Lateral  Distribution  of 
Longitudinal  Perturbation  Velocity 
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Figure  17.  Effect  of  on  the  Lateral  Distribution 

of  Longitudinal  Perturbation  Velocity 
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Figure  18 


Depth  Averaged  Longitudinal  Velocity 
Perturbations  at  y/b  =  0.9 


Depth  averaged  vertical  velocity  at  y/b=1.0 
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Figure  19. 


Depth  Averaged  Vertical  Velocities 
at  y/b  =  1.0 
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2.7  DEVELOPING  FLOW 
2.7.1  Formulation 

To  predict  longitudinal  variations  of  the  perturbation 
velocity  components,  the  full  set  of  Equations  163  to  163 
should  be  considered.  The  most  important  difference  from 
developed  flow  is  the  inclusion  of  terms  involving 
derivatives  with  respect  to  the  x  direction.  The  solution 
of  this  set  of  linear  equations  would  be  a  formidable  but 
not  impossible  task. 

A  simplified  solution  procedure  is  suggested  by  the 
following.  The  development  of  the  superelevation  and 
associated  inertial  effects  is  very  rapid  ( L  -  b)  and  is 
centered  about  a  change  of  curvature  (see  Appendix  A). 
Secondary  circulation  development  occurs  somewhat  more 
slowly  (Rosovskii,  1961)  but  is  still  very  rapid  compared  to 
the  shifting  of  the  longitudinal  velocity  profiles.  The 
behaviour  of  the  first  two  aspects  of  developing  flow  is 
also  somewhat  better  understood.  While  all  three  interact 
and  affect  each  other,  it  appears  that  the  first  two  reach 
developed  values  very  quickly  and  are  not  strongly  affected 
by  the  developing  longitudinal  velocity  distribution. 

This  analysis  of  developing  flow,  then,  considers  only 
Equation  163,  the  equation  for  the  vertical  vorticity  with 
the  following  simplifications: 

1.  The  terms  involving  v^  and  w^  are  to  be  evaluated 
using  the  fully  developed  solution  with  the  short 
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developing  region  represented  by  a  linear 
variation  of  the  term  from  zero  to  fully  developed 
between  x  =  0  and  x  =  0.1Cf/B. 

2.  The  term  involving  the  rate  of  change  of  curvature 

(1/R)  is  to  be  represented  as  a  triangular 

function  about  x  =  0  with  a  maximum  value  of  one 
and  a  base  extending  from  x  =  -1  to  x  =  1. 

3.  The  Ug  multiplying  the  first  term  and  the  rate  of 

change  of  curvature  term  are  taken  as  1.0  for 

simplicity. 

4.  The  ' vorticity'  n  is  taken  as 

3u-^ 

n  “  "  W~ 

neglecting  the  variation  of  v^  with  x  as  being 

small . 

5.  The  second  derivative  of  n  with  respect  to  x  in 

is  neglected  compared  to  the  vertical 
derivative.  This  simplification  reduces  the 

equation  to  parabolic  form  which  will  be  easier  to 
solve. 

The  equation  for  developing  longitudinal  velocity  is  then: 


ii  =  ill  rl5  +  ilni  +  v 

3x  1  „2  _2  0  J 
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The  boundary  conditions  in  y  and  z  are  the  same  as  for 
developed  flow  except  that  n  is  taken  simply  as  zero  at 
y  -  1  to  avoid  having  to  iterate  for  a  solution.  The 
velocity  u-^  is  then  derived  by  integrating  ri  with  respect  to 

y  • 


2.7.2  Solution 

Equation  204  may  be  solved  numerically  by  an  explicit 

technique  very  similar  to  the  successive-over-relaxation 

technique  used  to  solve  for  the  developed  flow  distribution. 
In  fact,  since  the  developed  values  of  v-j_  and  w^  are 

required,  it  is  only  necessary  to  modify  slightly  the 

program  used  for  developed  flow  and  to  save  the  longitudinal 

velocity  distribution  at  periodic  distance  steps. 

This  method,  however,  limits  the  maximum  distance  step 

for  a  stable  solution.  The  (conservative)  distance  step 

used  was: 


AX  =  0.25  — —  (Ay)2  (205) 

6vt 

where  Az  =  Ay. 

2.7.3  Comparison  to  Experiment 

This  problem  was  solved  for  parameters  corresponding  to 
Run  1  of  the  experimental  study  and  the  results  are  compared 
in  Figure  20.  This  plot  shows  the  variation  of  the  depth 


averaged  longitudinal  perturbation  velocity  at  y  =  0.9  with 
experimental  values  at  y  =  0.8  where  the  peak  occurred.  The 
graph  shows  a  negative  dip  at  the  beginning  of  the  bend  and 
then  an  exponential  climb  toward  the  fully  developed  value 
and  a  jump  at  the  end  of  the  bend  followed  by  an  exponential 
decay  toward  zero.  The  agreement  with  measured  values  is 
seen  to  be  reasonable. 

2.8  CONCLUSIONS 

A  theoretical  model  based  on  a  perturbation  of  the 
Reynold's  equations  for  a  moderate  curvature  ratio  is 
presented.  The  first  set  of  perturbation  equations  (linear 
in  curvature  ratio)  appear  to  predict  many  aspects  observed 
in  curved  channel  flows.  In  particular  it  is  possible  to 
calculate  an  approximate  redistribution  of  longitudinal 
velocity.  The  equations  derived  are  complex  but  linear  and 
are  solved  by  a  standard  numerical  technique  for  both 
developed  and  developing  flow.  Comparison  with  experiment 
indicates  reasonable  agreement. 

The  vertical  velocity  appears  to  have  an  important 
effect  on  the  redistribution  of  longitudinal  velocity.  The 
lateral  velocity  appears  only  to  smooth  the  distribution  out 
and  transfer  what  is  essentially  a  wall  effect  towards  the 
center  of  the  channel.  A  constant,  but  anisotropic,  eddy 
viscosity  was  used  to  model  the  turbulent  momentum  transport 
which  had  a  similar  effect.  It  is  felt  that  higher  order 
approximations  if  calculated  would  further  smooth  out  the 
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Figure  20.  Developing  Flow-Comparison  of  Theory 

to  Experiment 


longitudinal  velocity  distribution. 

The  weakest  part  of  the  analysis  is  likely  the  assumed 
form  for  the  zero  order  velocity  distribution.  While  a 
reasonable  approximation  over  the  upper  half  of  the  channel, 
it  seriously  underestimates  velocity  gradients  near  the  bed 
of  the  channel.  Other  distributions,  such  as  power  law, 
should  also  be  tried. 

The  analysis  is  restricted  to  rectangular  channels  with 
low  Froude  Number  flows.  The  second  restriction  is  not  too 
severe  but  for  application  to  practical  problems  the  problem 
should  be  generalized  to  more  arbitrary  geometries.  Perhaps 
the  easiest  method  would  be  to  allow  d  and  thus  3  to  be 
functions  of  x  and  y  and  perform  the  numerical  calculations 
using  nonconstant  coefficients.  This  would  also  result  in  a 
number  of  additional  terms  in  the  vorticity  equations  which 
would  have  to  be  accounted  for. 

In  summary,  the  analysis  presented  here  appears  to 
predict  many  features  of  curved  channel  flows.  It  needs 
more  verification  and  an  extension  to  more  arbitrary 
geometries  to  be  truly  useful  for  practical  problems. 
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3.  PART  THREE  -  EXPERIMENTAL  STUDY 

3.1  INTRODUCTION 

An  experimental  study  was  performed  to  provide  physical 
insights  into  the  problem  of  flow  in  curved  channels, 
validation  for  the  assumptions  and  results  of  the 
theoretical  analysis  and  to  add  significantly  to  the  body  of 
existing  experimental  data.  The  Laser  Doppler  Anaemometer 
(  LDA)  provides  unprecedented  accuracy  and  resolution  in 
liquid  flow  measurements  and  the  experimental  facility  was 
designed  specifically  for  the  DISA  LDA  system  in  the 
University  of  Alberta's  T.  Blench  Graduate  Hydraulics 
Laboratory. 

The  objectives  of  the  experimental  study  were  as 
follows : 

1.  The  physical  dimensions  of  the  flume  should  be  of 

realistic  proportions  and  the  channel  should  be  as 
large  as  space  permits. 

2.  The  total  angle  turned  by  the  bend  should  be  as 

large  as  possible  so  that  the  full  extent  of 

developing  flow  would  be  observed,  if  possible. 

3.  The  experiment  should  be  designed  to  take  full 

advantage  of  the  LDA's  unique  capabilities  for 
velocity  and  turbulence  measurements  including: 

(a)  detailed  longitudinal  velocity  measurements 
which  would  also  provide  longitudinal  bed 
stresses 
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(b)  detailed  lateral  mean  velocity  measurements 

which  would  be  measured  directly  and  which  in 
combination  with  the  longitudinal  velocity 

measurements  near  the  bed  would  provide  the 

angle  of  the  shear  stress  vector  at  the  bed. 

(c)  detailed  measurements  of  longitudinal  and 
lateral  turbulence  intensities. 

Due  to  the  detail  of  each  experiment,  only  two 
experiments  were  performed  and  only  the  depth  of  flow  and 

the  discharge  were  different  between  the  two  runs. 
Measurements  were  made  at  over  5,600  locations,  each 
involving  10,000  individual  velocity  resolutions  in  each 
direction,  for  a  total  of  over  one  hundred  million  velocity 
realizations . 

3.2  FACILITIES  AND  PROCEDURES 
3.2.1  Flume 

The  layout  of  the  flume  designed  and  used  for  these 
experiments  is  shown  in  Figure  21.  Also  included  in  the 
figure  are  the  locations  of  the  measurement  stations.  It  is 
comprised  of  a  straight  section  13.4  m  long  and  a  single 
bend  of  3.66  m  centreline  radius  over  an  arc  of  270° 
(distance  of  17.2  m)  followed  by  a  2.4  m  straight  exit 
section.  The  channel  was  1.07  m  wide  and  0.20  m  deep.  The 
slope  was  adjustable  but  was  set  at  1/1,200  for  both 
experiments.  At  this  slope  (necessarily  small  to  allow  a 
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Flume  Layout  and  Test  Locations 
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reasonable  Froude  number),  there  was  only  a  few  centimetres 
difference  between  the  bed  elevation  at  the  point  where  the 
channel  intersected  itself  so  the  flume  terminated  there  in 
a  drop  box. 

The  flume  was  constructed  of  galvanized  sheet  metal 
with  plexiglass  inserts  at  the  measurement  locations.  The 
bed  of  the  flume  was  a  minimum  of  1.5  m  above  the  floor  of 
the  lab  to  allow  access  from  underneath  by  the  LDA  system. 
The  supports  were  also  arranged  to  minimize  interference 
with  the  LDA  system. 

The  water  was  circulated  by  pump  from  a  level  regulated 
sump  to  the  bottom  of  a  head  box  equipped  with  an  overflow 
back  to  the  sump.  The  height  of  the  flume  allowed  the  sump 
and  head  box  to  be  above  floor  level. 

3.2.2  Laser  Doppler  Anaemometer 

The  LDA  technique  is  a  comparatively  recently  developed 
method  for  non-intrus ive ,  instantaneous  fluid  velocity 
measurements.  A  review  of  the  history,  theory  and  range  of 
application  of  this  technique  is  given  by  Durst,  Melling  and 
Whitelaw  (1978).  Buchave,  George  and  Lumley  (1979)  have 
provided  an  analysis  of  accuracies  to  be  expected.  The 
details  of  the  particular  2  colour  laser  equipment  used  is 
given  in  the  DISA  Manuals. 

The  laser  beams  entered  the  flow  from  the  bed  of  the 
flume  which  allowed  the  measurement  of  longitudinal  and 
lateral  velocity  components. 
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Figure  22  is  a  schematic  diagram  showing  the  main 
components  of  the  LDA  measurement  system  as  well  as  the 
control,  data  aquisition  and  analysis  systems.  The  laser 
was  a  4W  Argon- Ion  laser  with  special  optics  to  enhance  the 
blue  (488  nm)  line  and  an  etalon  to  improve  coherence. 
Actual  power  used  was  about  500  mW. 

The  LDA  optics  performed  the  dual  function  of  laser 
beam  manipulation  and  scattered  light  collection.  The 
single  beam  from  the  laser  was  split  into  three:  one  blue 

(488  nm) ,  one  green  (514.5  nm)  and  one  mixed  with  an 
approximate  power  distribution  of  25%/25%/50%.  The  two 
single  colour  beams  were  frequency  shifted  40  MHz  by  means 

of  a  Bragg  cell  to  enhance  signal  quality  and  insure  flow 

direction  information.  To  further  enhance  signal  quality 
and  to  reduce  the  size  of  the  measuring  volume/  a  beam 

expander  (1.9:1)  was  used. 

The  three  beams  exited  the  optics  parallel  to  and 
concentric  about  the  optical  axis.  The  optics  were  rotated 
so  that  the  blue  and  mixed  beams  were  perpendicular  to  the 
longitudinal  direction  while  the  green  and  mixed  beams  were 
parallel  to  the  longitudinal  direction. 

The  traversing  system  consisted  of  a  plane  mirror  and 
the  focusing  lens  mounted  on  threaded  tracks.  The  position 
of  the  mirror  and  lens  was  controlled  by  computer.  The 

nominal  position  accuracy  was  ±0.01  mm.  Focal  lengths  of 
the  lenses  used  were  80  mm  and  160  mm.  A  benefit  of  this 
system  was  that  the  laser  and  optics  remained  stationary 
during  traversing. 
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Figure  22 


Schematic  Diagram  oE  LDA  Setup 
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The  traversing  mechanisms  were  limited  to  0.6m  total 
travel  which  was  more  than  adequate  in  the  vertical 
direction  but  covered  only  about  half  the  width  of  the 
channel.  The  entire  LDA  system  was  mounted  on  a  trolley 
which  was  manually  positioned  twice  at  each  measurement 
section.  A  string  through  the  centre  of  channel  curvature 
was  used  to  align  the  traverse  in  a  radial  direction. 

The  flow  was  seeded  with  latex  paint  in  concentrations 
of  less  than  2  ppm.  This  material  was  found  to  be  very 
effective  at  a  minimum  cost.  Seeding  was  essential  for 
measurements  with  the  very  short  focal  length  lenses  where 
the  measuring  volume  was  very  small. 

The  scattered  light  was  collected  by  the  same  optics 
and  separated  by  colour  sensitive  mirrors.  The  light  was 
converted  to  electrical  signals  by  photomultipliers  and 
passed  through  the  dual  channel  channel  frequency  shifter  to 
two  identical  counters.  The  counters  evaluated  the  Doppler 
frequency  for  each  burst  and  this  value  was  transmitted  to  a 
computer  for  analysis.  A  coincidence  filter  was  used  to 
assure  approximately  simultaneous  information  on  each 
channel.  One  of  the  counters  also  provided  the  total  number 
of  fringes  in  the  burst.  This  assured  only  one  reading  per 
burst  and  was  also  useful  in  the  analysis  of  the  data. 

3.2.3  Data  Reduction 

The  coincident  pairs  of  Doppler  frequencies  were  stored 
as  they  arrived  in  the  computer  memory.  Due  to  the  limited 
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memory  capacity,  only  2,000  pairs  could  be  stored  at  one 
time.  These  were  then  reduced  and  successive  sets  of  data 
taken  until  the  desired  sample  size  was  obtained.  For  these 
experiments  a  sample  size  of  10,000  was  used.  The  average 

data  rate  was  about  200  samples  per  second. 

The  flow  statistics  obtained  were:  un ,  the  mean 

velocity  measured  by  the  green  set  of  beams;  un ,  the  mean 

D 

“  2 

velocity  m  the  blue  direction;  u*  ,  the  mean  square 

Cj  _ 

2 

turbulence  intensity  in  the  green  direction;  u'  ,  the 

D 

mean  square  turbulence  intensity  in  the  blue  direction;  and 
u'  u*  ,  the  covariance  of  the  velocity  fluctuations. 

These  values  were  calculated  for  Run  1  using  a 

two-dimensional  bias  correction  method  (Buchave,  George  & 

Lumley,  1979).  The  formulas  used  were: 
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uB,  =  velocity  reading  from  blue  beams 
Uq„  =  velocity  reading  from  green  beams 


N  =  number  of  accepted  samples 


(211) 


The  data  from  the  second  run  was  analyzed  using  the 
total  number  of  fringes  in  a  burst  which  is  related  to  the 
particle  transit  time  by: 


(211) 
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where 


tp  =  particle  transit  time 

Nf  =  number  of  fringes  observed  (zero  crossings) 

i 

fD  =  doppler  frequency 

The  averages  were  then  computed  using  tp  as  a  weighting 
coefficient.  For  example: 


UG 


N 

Z 

i  =  l 


uGi  Pi 


N 

I 

i  =  l 


(212) 


P. 

1 


The  two  methods  gave  almost  identical  results  for  the  same 
profiles  but  the  second  method  allowed  much  faster 
computation. 

The  procedure  used  a  two  pass  method  to  eliminate 
extreme  outlying  values.  Preliminary  statistics  were 
calculated  using  a  sample  size  of  300  points.  The 
calculation  then  proceeded  on  the  full  sample  with  any 
points  more  than  4  standard  deviations  from  the  mean  being 
excluded.  If  more  than  2.5%  of  the  points  were  thus 

excluded  the  calculation  was  repeated  once  more  using  the 

new  mean  and  standard  deviation.  Generally,  less  than  1.0% 

of  the  readings  were  excluded. 

It  was  also  necessary  to  correct  the  vertical 

coordinate  for  refraction.  The  following  formula  was  used: 


z  =  (  z  '  -  z  '  T ) 
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where : 


z 


z  1 

z  1 


I 


n 


distance  from  the  inside  of  the  bed  to 
measuring  point, 
traverse  reading. 

is  traverse  reading  when  measuring  point  is 
focused  on  the  inside  of  the  wall, 
refractive  index  of  water  (1.33) 
beam  intersection  angle 


For  each  location  the  coordinates  and  the  5  statistics 
were  recorded  on  floppy  disk  for  later  retrieval  and 
analysis. 

The  data  was  taken  in  the  form  of  vertical  profiles 
starting  as  close  to  the  bed  as  possible  and  proceeding 
upward  to  the  surface.  Eleven  profiles  were  taken  at  each 
section  with  about  twenty  points  on  each  profile.  The 
spacing  of  the  points  was  about  0.1  mm  near  the  bed  to  about 
6-10  mm  through  the  upper  part  of  the  flow.  The  approximate 
bed  location  was  noted  at  the  start  of  each  profile.  It  was 
not  possible  to  make  measurements  near  the  water  surface 
because  of  reflections  from  the  undulating  free  surface. 

3.3  EXPERIMENTAL  RESULTS 

The  significant  details  of  the  two  experiments 
performed  are  given  in  Table  4. 


' 


Table  4.  Significant  details  of  the  experiments 
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The  entire  set  of  data  for  Run  1  is  reproduced  in 
Figures  23  to  87.  The  information  of  u  velocity,  v 
velocity,  u'  fluctuation  intensity,  v'  fluctuation 
intensity  for  each  of  thirteen  sections  are  shown  on 
separate  plots.  An  additional  plot  for  each  section  shows 
the  u  velocity  with  the  vertical  scale  magnified  focussing 
on  the  region  very  close  to  the  bed.  All  dimensions  are  in 
millimetres. 

The  figures  are  organized  to  show  the  variation  of  the 
quantity  around  the  bend  so  the  u  profiles  are  shown  in 
sequence  followed  by  expanded  u  profiles,  lateral  velocity 
profiles  and  turbulence  profiles.  The  second  run  is  not 
shown  since  the  parameter  differences  from  the  first  run  are 
not  large  enough  to  result  in  any  significant  qualitative 
difference.  All  of  the  data  for  both  runs  is  available  from 
the  writer  in  the  form  of  tables  or  on  computer  tape. 

The  u  velocity  profiles  in  particular  show  an 
interesting  variation.  The  first  two  sections,  before  the 
start  of  the  bend,  have  an  almost  uniform  lateral 
distribution  with  only  the  two  outermost  profiles  being 
noticeably  different.  The  next  profile,  just  past  the  start 
of  the  bend  and  the  one  after  it  show  the  velocity  larger  on 
the  inside  of  the  channel  than  on  the  outside.  The  furthest 
outside  profile  is  already  catching  up  to  the  rest  by  the 
fourth  section,  however.  The  fifth  and  sixth  sections  show 
the  two  inside  profiles  falling  off  while  the  rest  become 
almost  uniform  again  (the  sixth  section  is  90°  into  the 
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bend).  Subsequent  profiles,  up  to  the  eleventh  show  a 
consistent  spreading  out  of  the  profiles  with  the  velocity 
increasing  at  the  outside  and  decreasing  at  the  inside  of 
the  bend.  The  twelfth  section,  just  past  270°  (the  end  of 
the  bend)  shows  a  sudden  increase  of  velocity  skewness  from 
the  eleventh  section  and  the  thirteenth  shows  an  even  larger 
variation.  The  shape  of  the  velocity  profiles  in  the  latter 
sections  of  the  bend  are  also  significantly  different  from 
the  corresponding  straight  channel  profiles.  They  are 
almost  uniform  and  in  some  cases  receding  through  the 
greater  part  of  the  channel  whereas  in  the  straight  channel 
the  profiles  increase  throughout  the  depth. 

The  expanded  u  velocity  profiles  indicate  much  the  same 
variation  over  the  bend  as  the  full  profiles,  although  the 
slopes  of  the  velocity  profiles  near  the  bed  may  be  taken  as 
indicative  of  the  bed  stress  at  that  location.  Many  of 
these  profiles  appear  not  to  pass  through  the  origin, 
indicating  a  small  error  in  fixing  the  bed  location.  This 
error  is  less  than  a  few  tenths  of  a  millimetre. 

The  lateral  velocity  (v)  profiles  show  the  development 
of  the  secondary  circulation  around  the  bend.  The  profiles 
in  the  straight  channel,  especially  the  second  one,  also 
show  evidence  of  the  much  smaller  secondary  flow  expected  in 
a  straight  channel.  The  maximum  magnitude,  near  the  surface 
at  the  furthest  outside  profiles,  directed  toward  the  centre 
of  the  channel  and  approximately  symmetrical  is  about  two 
per  cent  of  the  channel  mean  velocity. 


. 


At  the  third  section,  relatively  large,  inward 
corrected,  lateral  velocities  are  observed.  This 

corresponds  to  the  shifting  of  the  longitudinal  velocity 
maximum  towards  the  inside  here.  At  the  fourth  section  the 
circulation  profile  is  much  more  developed  but  now  there  is 
clearly  a  large  mean  outward  component.  Again  this 

corresponds  to  a  shifting  longitudinal  velocity  profile.  By 
the  fifth  section,  the  secondary  circulation  is  fully 
developed  and  in  fact  appears  to  actually  decline  in 

magnitude  through  the  rest  of  the  bend.  The  furthest 
outside  profile  also  shows  evidence  of  a  counterrotating 

circulation  cell.  The  twelfth  section  also  shows  a  net 

outward  discharge  and  the  thirteenth  section,  in  the 
straight  reach  following  the  end  of  the  bend  shows  the 
secondary  circulation  decaying. 

2 

The  turbulence  intensity  profiles,  first  u'  and  then 

2 

v'  are  also  presented  although  interpretation  is  not 
simple.  The  measurements  of  u'w'  turned  out  to  be  too 
small  to  show  any  systematic  variation  at  all.  Since  this 
correlation  corresponds  to  the  turbulent  stress  that  would 

act  on  the  side  walls  it  was  expected  to  be  very  small 

through  the  main  central  part  of  the  channel.  In  fact,  it 
was  smaller  than  the  experimental  scatter. 
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Figure  40.  Expanded  u  Velocity  Distribution 


146 


o  cd  o  on  co 


CO  CM  CM  'OOtDCOOCDO 
^  CO  CsJ  — '  •  O  «-<  CsJ  CM  OO 
!  I  1  I  i  O  — <  CM  CO  -3-  '<3- 

II  II  II  U  II  II  II  M  II  II  II 


B04+X*4-KN>»X 


y  WKe  a 

K 


****«■© 

K 


e 

K 


»««e 

K 


-  CO 


-  LO 


K 


hMI-KI 


X  >-  X 
K 


1 - 1 — 

oe-  o 


— i - 1 - 1 - i - r 

02-o  oro 

(  s/uj  )  n 


co 

LiJ 

^  cn 

h~ 

LU 

n 

H 


-  cn 


-  CM 


H 


N 


0 


— i — 
0*  •  0 


00*0 


Figure  41.  Expanded  u  Velocity  Distribution 
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Figure  46.  Expanded  u  Velocity  Distribution 
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Figure  54.  v  Velocity  Distribution  (x 
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Figure  68.  u'  Turbulence  Intensity  Distribution 
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Figure  79.  v'^  Turbulence  Intensity  Distribution 

(x  =  3.83  m) 
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Figure  83.  v*  Turbulence  Intensity  Distribution 
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Figure  84.  v’  Turbulence  Intensity  Distribution 
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Figure  36.  v'  Turbulence  Intensity  Distribution 
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Figure  87.  v*^  Turbulence  Intensity  Distribution 
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3.4  Analysis  of  Results 

The  results  of  the  two  experiments  were  reduced  to 
profile  scales  to  aid  in  evaluation  of  the  overall  trends  of 
the  developing  flow.  The  mean  longitudinal  and  lateral 
velocity  for  each  profile  were  computed  by  direct 

integration.  The  variation  of  these  quantities  would  show 
the  redistribution  of  velocity  to  the  inside  of  the  bend  at 
the  beginning  and  then  towards  the  outside  through  the  rest 
of  the  bend. 

The  longitudinal  bed  stress  was  derived  by  analysis  of 
the  u  velocity  profile.  A  first  estimate  of  u*  was  obtained 
by  fitting  a  semi-log  line  through  the  bottom  10  points  of 
the  profile.  This  estimate  was  then  used  to  eliminate  any 
points  less  than  70  v/u*  from  the  bed.  The  regression  was 
done  again  with  the  refined  set  of  points  until  there  was  no 
change  in  u*. 

The  total  bed  stress  vector  was  assumed  to  be  in  the 
same  direction  as  the  velocity  vector  near  the  bed  ( Nash  and 
Patel,  1972).  Polar  plots  of  velocity  (Figure  88,  for 
example)  indicated  that  in  many  profiles  there  were  several 
points  measured  for  which  this  appeared  to  be  true.  The  bed 
stress  angle  was  derived  by  a  linear  regression  of  these 
points.  The  criteria  used  for  the  upper  limit  was  that: 

<  12.0  (214) 

u* 


as  suggested  by  Nash  and  Patel  (1972).  Some  profiles  in  the 


✓ 


' 


v  (m/s) 

.06  -0.04  -0.02  0.00  0.02  0.04  0.06 


194 


a 

© 

▲ 

+ 

x 

♦ 

♦ 

X 

z 

Y 

X 


+A 


Y/B 

-0.45 

-0.40 

-0.30 

-0.20 

-0.10 

0.00 

0.10 

0.20 

0.30 

0.40 

0.45 


X  x 

X 

X 


♦ 

i 


a  ▲ 
▲ 


+ 

+ 

+ 


* 

♦  x 

X 


®_ 

o 


^JG! 


Z 

Z 


z 

*  z  ^ 

X 

z  ** 

*  X 


ex 


H  a 


+®° 


3?  *  D  .  _  ® 

♦  +  j?®  ® 

X  Xa 


IS 


E 


x  r 
x 
z 


X 

X 


V 

X  z 


X  ♦ 


x  A 

X  .  ..  + 


+  i  x+ 


Y 

z  * 


♦  X  * 


x  XX, 
X  * 


♦ 

Y  Z 
Y 


*  ♦ 


Y 


♦ 

Z 


+  ♦  y*.  Y+  ^ 


X  X 


z  z  X 
X 


O 

I 


oo 

o 


o 

I _ 

0.00 


0.10 


— i - 1 - 1 - 1 - 1 - 1 - 

0.20  0.30  0.40  0.50 

u  (m/s) 


Figure  88 


Polar  Plot  of  Velocity  Vector 
( x  =  11.49  m) 


X  X 


second  run  contained  no  points  in  this  region  and  no  bed 
stress  angle  was  computed. 

Once  the  bed  stress  angle  was  obtained  it  was  a  simple 

matter  to  calculate  the  lateral  stress  magnitude  as: 

; 

Tv0  2 

=  tan<*>  u*  (215) 


To  obtain  a  measure  of  the  secondary  (circulatory)  flow 
the  following  integral  was  also  calculated: 

__  1  d 

vz  =  ^  /  z(v  -  v)  dz  (216) 

This  moment  of  the  lateral  velocity  distribution  gives  an 
indiction  of  the  relative  strength  of  (zero  depth  average) 
circulation. 

All  of  these  quantities  were  calculated  by  a  computer 
program  which  is  included  in  Appendix  B  for  reference.  The 
regression  curve  fitting  subroutine  used  was  from  Peterson 
and  Howells  (1981). 

3.5  DISCUSSION 

Figures  89  and  90  show  the  longitudinal  variation  of 
depth  averaged  longitudinal  velocity  for  the  two  runs.  From 
a  relatively  uniform  lateral  variation  before  the  beginning 
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Figure  89.  Longitudinal  Variation  of  Depth  Averaged 

Longitudinal  Velocity  -  Run  1 
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Figure  90.  Longitudinal  Variation  of  Depth  Averaged 

Longitudinal  Velocity  -  Run  2 
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of  the  bend  (x/b  -  5)  the  distribution  skews  inward  for 
the  first  part  of  the  bend.  After  x/b  of  about  10  the 
distribution  starts  to  skew  outward  and  this  trend  appears 
to  continue  for  the  rest  of  the  bend.  At  the  end  of  the 
bend  (x/b  =  32)  a  further  sudden  outward  skewing  takes 
place.  Due  to  the  short  exit  channel  available,  no 
significant  return  to  a  distribution  was  observed. 

The  second  run  appears  to  be  even  less  developed  by  the 
end  of  the  bend  than  the  first  run.  The  furthest  inside 
profile  (at  b/y  =  -0.8)  in  both  runs  shows  a  larger  velocity 
decrease  than  the  corresponding  outside  profile  (y/b  =  0.8) 
shows  a  velocity  increase.  Both  runs  also  indicate  a  slight 
decrease  in  cross-sectional  average  velocity  through  the 
bend  but  a  small  rebound  after  the  bend.  The  ninth  and 
tenth  section  appeared  to  behave  slightly  anomalously,  with 
lateral  variation  appearing  too  large  at  section  nine  and 
too  small  at  section  ten.  An  uneven  bed  in  this  vicinity 
may  be  the  explanation. 

Figures  91  and  92  show  the  longitudinal  variation  of 
the  longitudinal  bed  stress  in  the  form  of  a  shear 
velocity.  With  somewhat  more  scatter,  possibly  reflecting 
the  calculation  method  more  than  anything  else,  these 
distributions  basicly  follow  the  pattern  of  the  depth 
averaged  longitudinal  velocity.  The  magnitude  of  the  cross 
channel  variation  is  relatively  larger.  Sections  nine  and 
ten  again  appear  somewhat  out  of  line  with  the  rest  of  the 


sections . 


Figures  93  and  94  show  the  variation  of  bed  stress 
angle  as  derived  from  the  polar  plots.  These  distributions 
show  a  significant  variation  across  the  channel  with  angles 
at  both  y/b  =  ±0.8  about  60%  of  the  more  central  profiles. 
The  distributions  appear  to  be  approximately  constant  after 
x/b  -  5  decreasing  slightly  in  magnitude,  if  anything.  The 
maximum  angle  for  Run  1  ( tan<j>  =  0.20)  corresponds  to  a  value 
of  K  of  about  12  with  an  average  value  of  about  10  in  the 
central  part  of  the  channel  through  most  of  the  bend. 
Insufficient  profiles  meeting  the  criteria  were  available  to 
make  any  comments  for  run  2. 

Figures  95  and  96  show  the  magnitude  of  the  lateral  bed 
stress  as  calculated  from  Equation  215.  These  variations 
are  similar  in  general  form  to  the  bed  stress  angle  plots 
but  appear  to  vary  more  smoothly.  The  lateral  variation  is 
even  more  pronounced  however  and  the  lateral  stress  at 
y/b  =  0.4  is  clearly  the  largest  for  the  greatest  part  of 
the  bend. 

Figures  97  and  98  show  the  variation  of  the  depth 
averaged  lateral  velocity  along  the  channel.  These  values 
are  generally  small  (less  than  5  per  cent  of  the  mean 
longitudinal  velocity)  but  show  an  interesting  variation. 
At  the  beginning  of  the  bend  (x/b  =  0)  there  is  a  relatively 
large  inward  flow,  followed  almost  immediately  by  a 
relatively  large  outward  flow.  For  Run  1,  the  flow  appears 
outward  for  most  of  the  channel  while  for  Run  2  the  flow 
appears  inward.  Since  the  longitudinal  velocity 
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Figure  91.  Longitudinal  Variation  of  Longitudinal 

Shear  Velocity  -  Run  1 
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Figure  92.  Longitudinal  Variation  of  Longitudinal 

Shear  Velocity  -  Run  2 


0.05 


202 


O 

o 

6 


in 

o 

d 


o 

d 


in 

o 

q 

O 

o 

<P  UD| 


Figure  93.  Longitudinal  Variation  of  Bed  Stress 

Angle  -  Run  1 
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Figure  95.  Longitudinal  Variation  of  Lateral  Bed 

Stress  -  Run  1 
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Figure  96.  Longitudinal  Variation  of  Lateral  Bed 
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Figure  97.  Longitudinal  Variation  of  Depth  Averaged 

Lateral  Velocity  -  Run  1 
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Figure  98.  Longitudinal  Variation  of  Depth  Average 

Lateral  Velocity  -  Run  2 


distribution  continued  to  skew  outwards  in  Run  2  (even  more 
than  in  Run  1)  it  may  be  concluded  that  a  small  alignment 

error  was  present  in  Run  2  causing  the  lateral  velocity 
readings  to  appear  more  negative  than  they  were. 

The  last  two  plots,  Figures  99  and  100  show  the 

variation  of  the  first  moment  of  the  lateral  velocity 
distribution.  Indicating  the  relative  strength  of  the 
circulation  at  this  point,  these  plots  show  the  circulation 
developing  very  quickly  (by  about  x/b  =  5).  They  also 
indicate  that  the  circulation  also  declines  noticeably  over 
the  length  of  the  bend.  This  is  more  apparent  in  Run  2. 
The  circulation  near  the  side  walls  is  considerably  less 
(about  half)  than  in  the  center  of  the  channel.  The 

beginning  of  the  decay  after  the  end  of  the  curve 

(y/b  =  32.0)  is  also  shown. 


3.6  CONCLUSIONS 

Two  complete  experimental  runs  were  performed.  The 
Laser  Doppler  Anaemometer  provided  velocity  measurements  in 
both  the  longitudinal  and  lateral  directions  of 
unprecedented  detail  and  precision.  The  time  involved  for 
each  run  precluded  a  wider  range  of  geometrical  flow 
parameters  however. 

Despite  a  total  included  angle  of  270° ,  neither  run 
exhibited  fully  developed  flow.  The  secondary  circulation 
developed  quickly  but  the  longitudinal  velocity  distribution 
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Figure  99.  Longitudinal  Variation  of  Secondary 

Circulation  Intensity  -  Run  1 
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Figure  100.  Longitudinal  Variation  of  Secondary 

Circulation  Intensity  -  Run  2 


took  far 

longer . 

The 

longitudinal  and  lateral  velocity  profiles  were 

detailed 

enough  near  the  bed  to  allow  a  direct  estimate  of 

the  velocity  vector  angle  near  the  bed.  This  provided  a 
reliable  estimate  of  the  bed  stress  angle  and  the 
controversial  constant  'K'.  It  was  found  that  this  angle 
and  constant  remained  approximately  constant  through  the 
bend  but  varied  considerably  across  the  channel.  A  typical 
value  of  K  near  the  centre  of  the  channel  was  about  10, 


which  is 

in  close  agreement  with  the  assumption  of  Part  One. 

The 

theoretical  analysis  of  Part  Two  suggests  that 

vertical 

velocities  should  be  measured,  especially  in  the 

vicinity 

of  the  outside  wall  and  also  the  lateral  turbulent 

shear  stress  in  this  vicinity.  The  first  should  be  measured 
because  of  it's  importance  in  the  skewing  of  the 
longitudinal  velocity  profile  and  second  because  of  it's 


relation 

to  the  lateral  eddy  viscosity  coefficient  which 

also  was 

important  in  this  prediction.  These  two  quantities 

are  very 

small  and  would  require  extreme  care  even  with  the 

best  facilities.  Only  a  few  sections  would  be  necessary  to 
validate  or  invalidate  the  theoretical  assumptions  made.  It 
would  also  be  desirable  to  measure  the  distributions  of 
u '  w '  and  v '  w '  directly  to  provide  better  stress 
information  that  does  not  rely  on  assumed  velocity 
distributions.  These  again  would  require  special  facilities 
and  care.  The  lateral  stress  v'w'  especially  would  be 


difficult  to  measure  because  of  geometrical  problems. 


Finally  it  would  be  desirable  to  make  detailed  measurements 
for  a  wider  range  of  geometrical  and  flow  parameters.  The 
bed  and  side  roughness  especially  could  be  changed,  although 
the  use  of  the  LDA  would  be  impaired  if  a  clear  viewport  was 
not  available. 

In  conclusion,  detailed  experiments  were  performed 
which  further  illustrate  the  flow  mechanics,  provide  a 
reference  for  comparison  with  theoretical  predictions  and 
add  significantly  to  the  body  of  experimental  knowledge. 
There  remains,  however,  several  apsects  of  problem 
completely  unexplored  and  a  wide  range  of  parameters 


unmeasured . 


CONC LUS ION 


In  this  study,  the  present  knowledge  was  reviewed  with 
the  objectives  of  applying  it  to  practical  problems  and 
identifying  directions  for  further  research.  The  lack  of 

prediction  methods  for  the  developing  and  developed 
longitudinal  velocity  distribution  was  identified  as  the 
most  important  gap  in  the  existing  knowledge.  It  was  also 
noted  that  few  detailed  experiments  had  been  performed  and 
many  important  quantities  such  as  lateral  bed  stress  angle 
and  vertical  velocity  components  and  turbulence  had  been 
measured  only  indirectly,  if  at  all. 

On  the  positive  side,  it  was  found  that  a  good 

understanding  of  many  aspects  of  curved  channel  flow  such  as 
superelevation  and  secondary  flow  could  be  applied  to 
practical  problems.  In  particular  it  is  possible  to 

estimate  the  maximum  scour  depth  in  an  alluvial  channel. 
The  method  is  still  subject  to  considerable  variation  but 
certainly  takes  better  account  of  the  physics  of  the  problem 
than  the  existing  regime  methods. 

In  an  attempt  to  predict  longitudinal  velocity 
variations  a  perturbation  analysis  was  performed  on  the 
governing  equations  of  flow.  Using  the  ratio  of  channel 
half  width  to  centreline  radius  of  curvature  as  the 
pertubation  parameter  it  was  found  that  a  first  order 

(linear)  approximation  included  most  of  the  interesting 
curved  channel  effects.  The  analysis  depended  on  an 
assumption  for  the  zero  order  (straight  channel)  velocity 
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distribution  and  this  may  be  the  greatest  weakness  in  this 
analysis. 

A  numerical  solution  for  fully  developed  flow  was  run 
for  a  range  of  channel  geometries  and  roughnesses.  This 
provided  an  estimate  of  the  maximum  bank  attack  velocity 
likely  to  be  encountered  in  a  given  problem.  A  simplified 
model  of  developing  flow  also  gave  reasonable  agreement  to 
experimental  data. 

The  method  appeared  to  account  for  all  the  velocity 
redistribution  factors  but  is  limited  to  rectangular 
channels.  A  useful  extension  of  this  work  would  be  to 
generalize  it  to  arbitrary  cross-sectional  shapes  and  also 
evaluate  different  initial  velocity  distribution 
assumptions.  Less  important  but  interesting  would  be  to 
investigate  different  turbulence  closure  assumptions. 
Calculation  of  a  second  order  approximation  would  also  be 
interesting  but  would  involve  considerable  effort. 

An  experimental  study  comprising  detailed  measurements 
of  longitudinal  and  lateral  velocity  components  for  two  flow 
situations  was  performed.  The  intricacies  of  the  flow  were 
revealed  by  these  experiments  and  the  results  were  used  to 
validate  the  theoretical  predictions.  The  long  bend  (270° 
included  angle)  covered  most  of  the  developing  flow  region 
but  developed  flow  was  not  observed.  Bend  entrance  and  exit 
effects  were  observed.  Measurements  of  vertical  velocities 
and  turbulent  shear  stresses  have  yet  to  be  made  and  it 
would  be  desirable  to  obtain  detailed  measurements  for  a 
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wider  range  of  geometrical  and  flow  parameters. 

In  summary,  a  detailed  literature  review,  theoretical 
analysis,  and  experimental  study  have  been  performed  on  the 
subject  of  turbulent  flow  in  curved  channels.  Although  much 
more  research  needs  to  be  done,  it  is  felt  that  the 
direction  taken  by  this  study  has  provided  some  valuable 
results  and  the  potential  for  further  study. 
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APPENDIX  A 


WATER  SURFACE  AT  A  CHANGE  OF  CHANNEL  CURVATURE 

Introduction 

The  problem  of  flow  in  curved  open  channels  is 

important  in  hydraulic  engineering  and  has  been  the  subject 
of  much  recent  research.  In  most  of  this  work  the  lateral 
pressure  gradient  has  been  assumed  to  be  an  algebraic 

function  of  the  local  channel  curvature.  Clearly,  however, 
in  cases  where  the  curvature  changes  abruptly,  the  water 
surface  adjusts  continuously.  Inspection  of  experimental 
water  surface  profiles  indicates  that  an  adjustment  length 
about  equal  to  the  channel  width  is  required. 

Calculations  of  developing  secondary  flow  may  be  quite 
strongly  influenced  by  the  water  surface  configuration. 

Rosovskii  (1961)  and  Nouh  and  Townsend  (1979)  give  a 

development  length  (to  90%)  for  secondary  flow  as: 

L  -  2.0  HC f  (Al) 

where  H  is  the  average  depth  of  flow  and  C^  is  the 
non-dimensional  Chezy  coefficient.  Therefore  many 

situations  will  arise  where  this  length  is  comparable  to  the 
development  length  for  the  superelevation. 

Studies  which  do  account  for  this  phenomenon  such  as 
the  numerical  models  of  Leschizner  and  Rodi  (1979)  and  de 
Vriend  (1981)  or  potential  flow  solutions  such  as  Kamiyama 
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(1966)  neglect  free  surface  effects  and  therefore 
theoretically  apply  only  to  situations  of  small  Froude 
number.  It  is  thus  of  interest  to  investigate  the  effect  of 
Froude  number  on  the  water  surface  configuration. 

i 

This  paper  presents  a  simple  model  for  the  water 
surface  configuration  in  the  vicinity  of  an  abrupt  change  in 
channel  curvature  taking  free  surface  effects  into 
account.  The  flow  is  assumed  frictionless  and  the  curvature 
is  assumed  mild.  Preliminary  comparisons  with  experiments 
are  also  presented. 


Analysis 

Referring  to  Figure  A1  we  may  write  the  depth  averaged 
Euler  and  continuity  equations  for  a  rectangular  channel  in 
a  locally  cylindrical  coordinate  system  as: 
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wherein  the  centreline  radius  of  curvature,  R,  may  vary 


with  x. 


22  2 


X 


I 


Figure  A1 .  Definition  Sketches 
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The  pressure  distribution  has  been  assumed  hydrostatic 
and  all  friction  terms  have  been  neglected.  The  latter 
assumption  is  based  on  the  traditional  transition  assumption 
of  a  short  length  and  on  the  fact  that  even  for  fully 
developed  flow  in  a  curved  channel,  frictional  effects  on 
the  superelevation  are  small  (Yen  and  Yen,  1971). 

It  is  useful  to  non-dimensionalize  Equations  2  to  4 
using  the  undisturbed  straight  channel  mean  values  as 
scales. 

Let: 

u=Uu*;v=Uv*;h=  Hh*  (A5-A7) 


Also  let: 

x  =  bx*;  y  =  by*;  R  =  RQR*  (A8-A10) 

where  Rq  is  the  minimum  value  of  R  expected. 

These  may  be  introduced  into  Equations  2  to  4  and  after 
some  simplification  (starred  quantities  are  now  represented 
by  their  unstarred  equivalents)  we  obtain: 
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wherein : 


a  = 


( A14  ) 


and : 
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( A15) 


For  all  but  the  sharpest  bends  a  is  small  compared  to 
one.  This  leads  naturally  to  a  solution  of  the  above 
problem  in  terms  of  a  power  series  in  a.  Thus: 


u(x,y,FR,a)+  uQ(x,y,FR)+  au  1  (  x  ,y  ,FR)  +  a2u 2 ( x , y , F  ) +e tc . 

(A16) 

Also : 
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by  binomial  expansion. 

Equations  16  and  17  may  be  substituted  into  11  through 
12  and  then  grouped  by  degree  of  a  to  be  solved 
separately.  The  solution  to  the  zero  order  group  is  simply: 
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( A13-A20) 


The  first  order  group  with  18-20  substituted  becomes: 
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These  equations  may  be  solved  for  h-^  to  give 
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Boundary  conditions  at  y  =  ±1  may  be  obtained  from 

Equation  22: 
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The  boundary  conditions  for  Equation  28  are  simply: 


w  =  0  at  y  =  ±1  and  x  =  ±  «  (A29) 


Because  of  the  simple  boundary  conditions  the  boundary  value 
problem  represented  by  Equations  28  and  29  will  be  the  one 
actual  solved  with  u-^  and  h-^  given  later  by  Equations  21  and 
23. 

For  the  situation  of  a  long  straight  reach  followed  by 

a  long  curved  reach  7-7 — 7  may  be  specified  as  follows: 

k  (,  x ; 


=  0  ;  -°°  <  x  <  Xq  ( A30 ) 

=  1;  Xq<  x  <  « 

where  Xq  is  the  point  of  curvature  change. 

It  should  be  noted  that  Equations  28-29  form  a  linear 

l 

system  and  a  solution  for  arbitrary  p"(iry  ma^  £°rme<3  by 
superpositions  of  the  solution  for  Equation  28. 

Equation  28  is  elliptic  for  FR  <  1  and  hyperbolic  for 

Fr>  1.  Standing  cross  channel  waves  for  supercritical  flow 
are  correctly  predicted  (to  a  linear  approximation).  In  the 
following  only  the  subcritical  case  will  be  considered. 
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The  derivative  of  the  unit  step  function  defined  in 
Equation  30  is  the  Dirac  delta  function  denoted  as  D  (x  - 
Xq ) .  Letting : 


2 

a  = 


(1  -  F*> 


Equation  28  becomes: 
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This  boundary  value  problem  may  be  solved  by  the  Fourier 
Transform  method.  Take  the  Fourier  Transform  of  Equation  32 
with  respect  to  (x  -  Xq) 
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where  V  is  the  Fourier  transform  of  v^.  Rearrange  Equation 
33: 
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The  boundary  conditions  become: 
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The  general  solution  of  34  is: 

V  =  Acosh(^)  +  Bsinh(^)  +  c 
Applying  the  boundary  conditions: 
0  =  Acosh(j)  -  Bsinh  (^)  +  c 
0  =  Acosh(-)  +  Bcosh(-)  +  C 


Solving : 


B  =  0 


A  = 
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C  =  - 
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Thus  : 
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Take  the  inverse  Fourier  Transform 
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The  integral  may  be  evaluated  using  contour  integration  in 
the  complex  plane.  The  formula  is: 

/  (  )  =  £  (Residues  at  singularities  (A44) 

c 

enclosed  within  C) 


The  contour  used  was  along  the  real  axis  from  -°°  to  «  then 
along  a  semicircular  arc  at  |k|  =  «  back  to  k  =  -«  (real). 
The  singularities  enclosed  are  simple  poles  at 
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The  residue  at  each  pole  is  thus: 
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This  solution  may  also  be  verified  by  direct  substitution. 
The  corresponding  solutions  for  u^  and  h^  are: 


ui  = 


l 

n=0 


4  ( -1 ) 


n 


2  2 

(2n+l)  tt 


sm 


( 


( 2n+l ) uy 


)  A 


(  A43) 
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CO 


I 

n  =  0 


4  (-l)n 

(  2n  +  l )  2tt 2 


sin 


(lip-tim)  A 


where : 


A 


exp  (-Until  7Ta(  Xq-x)  )  ;  x  <  xQ 
2-exp  (~2n^  7ra(xQ-x));  x  >  xQ 


(  A49) 


( A50 ) 


and  : 


a 


(  1-F 


-1/2 


(A51) 


A  characteristic  length  may  be  derived  by  considering 
the  behaviour  of  this  function  at  y  =  1.  Define  L  to  be 

•D 

the  distance  between  the  points  where  h^'  is  5%  and  95%  of 
it's  final  value.  Lg  may  then  be  solved  from: 


0.05  =  I 

n  =  0 


(  2n  +  l )  2  7t 2 


( —  (  2  n+ 1 )  (  s  ^  'i 

exp  [ - 5 - -  77  ( - JJ 

z  2b 


( A52 ) 


Solving  numerically: 


2  i/2 

Ls  =  2.66b  (1  -  Fz) 


( A53) 


EXPERIMENTS 

To  provide  experimental  verification  for  the  preceeding 
analysis  a  few  preliminary  experiments  were  performed. 
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Existing  experiments  were  found  unsuitable  due  to  lack  of 
detail  in  the  transition  region. 

The  first  experiment  was  done  in  a  flume  with  a  single 
270°  bend.  The  average  flow  depth  was  61  mm;  the  average 
velocity  was  0.38  m/s ;  the  width  was  1.07  m  and  the 
centerline  radius  was  3.66  m.  The  Froude  number  was 

therefore  equal  to  0.495  and  a  was  0.146. 

The  water  surface  elevations  were  measured  with  a 
screwdriver  probe  and  a  sloping  manometer.  Estimated 
accuracy  was  ±  0.2  mm. 

Figure  A2  shows  a  comparison  of  Equation  49  and  the 
experimental  results.  The  variable  h*  is  given  by: 


(h  -  H) 
(U2b/gRQ) 


(A54) 


The  agreement  is  seen  to  be  reasonable. 

The  second  experiment  was  performed  in  an  'S'  shaped 
flume  where  the  curvature  was  abruptly  reversed.  In  this 
case  F  =  0.50  and  a  =  0.140.  Figure  A3  shows  a  comparison 
of  this  experiment  with  theoretical  curves  obtained  by 
superposition.  Again  the  agreement  is  observed  to  be 
reasonably  good. 


CQNC  LUS  IONS 

An  approximate  relation  for  the  water  surface 
configuration  at  an  abrupt  change  of  curvature  has  been 
derived  for  a  cnannel  of  mild  curvature.  A  simple 


1.25 
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Curved  Channel 
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Figure  A3.  Water  Surface  at  a  Sudden  Reversal  of 

Channel  Curvature 


expression  for  the  length  of  such  a  transition  as  a  function 
of  Froude  number  has  also  been  presented.  The  predicted 
water  surface  topography  was  shown  to  agree  reasonably  well 
with  a  few  preliminary  experiments. 


APPENDIX  B 


COMPUTER  PROGRAMS 


PROGRAM  CRVCNL 

The  computer  program  CRVCNL  is  a  numerical 
implementation  of  the  Engelund  (1974)  method  of  flow  and  bed 
topography  calculation  for  a  curved  channel.  The  method  has 
been  generalized  to  allow  an  arbitrary  meander  geometry. 
The  channel  curvature  is  calculated  from  the  centreline 
co-ordinates  of  the  stream  plan  supplied  by  the  user.  The 
program  outputs  velocity  (longitudinal)  corrections  and  bed 
depths  below  the  water  surface  as  well  as  various 
intermediary  calculated  quantities. 

This  program  is  intended  for  the  demonstration  of  the 
Engelund  method  only.  Included  are  input  instructions  and  a 
listing  of  the  FORTRAN  source  code. 

INPUT  TO  CRVCNL 


FORTRAN 
Card  #1 

1. 


2. 

3. 


input  unit  =  5 

(5F10.4) 

Average  depth  of  flow  for  a  straight  reach 
(consistent  length  unit); 

Average  width  of  stream  (consistent  length  unit); 
Friction  factor  of  stream  defined  from: 


f 


(Bl) 
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4. 


Sediment  friction  angle  (in  degrees)  (  =  25  to  30 

degrees,  may  be  calibrated); 

5.  Sediment  transport  coefficient  (=4.0,  may  be 
calibrated) . 


Card  #2 
1. 

2. 

3. 


(  213 , F10 . 4 ) 

Number  of  longitudinal  sections  (NSP); 

Number  of  calculation  points  across  channel  (  NNP) ; 
Distance  between  sections  (consistent  length  unit); 


Card  #3  to  Card  #(NSP+2) 

1.  X  co-ordinate  of 
unit ) ; 

2.  Y  co-ordinate  of 
unit ) . 


( 2F10. 4) 

channel  plan 


channel  plan 


(consistent  length 


(consistent  length 


These  co-ordinates  are  determined  by  choosing  equally 
spaced  points  along  the  channel  and  referring  each  point  to 
an  arbitary  Cartesian  co-ordinate  system  with  the  positive 
x-direction  corresponding  roughly  to  downstream  along  the 
valley  axis.  All  lengths  output  will  be  expressed  in  the 
same  unit  used  in  the  data  input. 


OUTPUT 


FORTRAN  output  unit  =  6 

The  output  of  the  various  variables  is  similar.  The 
longitudinal  co-ordinate  is  printed  on  the  left  of  the  page 
while  the  cross  channel  variation  of  the  variable  is  printed 
across  the  page.  The  distance  between  these  points  is  the 
width  divided  by  (NNP-1). 

PROGRAM  EXAMPLE 

The  program  is  demonstrated  on  a  situation  which  is  an 
idealization  of  the  North  Saskatchewan  River  near  the 
Mayfair  Golf  Course  in  Edmonton.  The  river  is  represented 
as  a  single  180°  bend  with  a  centreline  radius  of  750  m. 
Flow  conditions  are  based  on  the  two  year  flood.  All  length 
dimensions  are  in  metres.  The  input  data  file  is  listed  on 
Figure  B1 ,  the  final  bed  topography  output  is  shown  on 
Figure  B2,  and  the  program  source  code  is  listed  on 
Figure  B3 . 

PROGRAMS  PERTRB  AND  DEVFLO 

These  programs  perform  the  numerical  analysis 
procedures  described  in  Part  Two.  PERTRB  does  the  developed 
flow  calculations  and  DEVFLO  does  the  developing  flow. 
These  programs  are  included  for  reference  only  and  are 
listed  in  Figures  B4  and  B5 


. 
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PROGRAM  AN Z DAT 

This  program  analyzes  the  data  as  described  in  part 
three.  It  is  listed  in  Figure  B6. 


239 


Figure  B1 . 


1 

3.4, 

170. ,0.0050,23 

2 

29 

7  131.0000 

3 

-655 

. , -750.  , 

4 

-524 

. , -750.  , 

5 

-393 

. , -750.  , 

6 

-262 

. ,-750.  , 

7 

-131 

. ,-750.  , 

8 

o.  ,- 

750.  , 

9 

130. 

,-739. , 

10 

257. 

,-705. , 

1  1 

375. 

,-650. , 

12 

482. 

,-575. , 

13 

575. 

,-482. , 

14 

650. 

, -375. , 

15 

705. 

,-257. , 

16 

739. 

,-130. , 

17 

750. 

,0.  , 

18 

739. 

, 130.  , 

19 

705. 

,257.  , 

20 

650. 

,375.  , 

21 

575. 

,482.  , 

22 

482. 

,575.  , 

23 

375. 

,650.  , 

24 

257. 

,705.  , 

25 

130. 

,739.  , 

26 

0. ,750.  , 

27 

-131 

. ,750.  , 

28 

-262 

. ,750.  , 

29 

-393 

. ,750. , 

30 

-524 

. ,750.  , 

31 

-655 

. ,750. , 

End  of  file 


CRVCNL  Example  Input 
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CURVED  CHANNEL  BED  TOPOCRAPMV  OUTPUT 

S(M)  BED  TOPOGRAPHY 


0 . 0 

©  .  34006*0 1 

© . 34ooe  *o i 

© . 3  4  O© E  *0 1 

© 

0  .  1  3 1 OE  +03 

©  . 34006*0 1 

0 . 34006*0 1 

0 . 3400E*0 1 

0 

0 . 26206*03 

0 . 34006*0 1 

© . 34006*0 1 

0 . 34006*0 1 

© 

0  39306*03 

© . 3  4  O©  E  *0 1 

0 . 34006*0 1 

©  34006*01 

0 

0 . 52406*03 

O . 46056*01 

O  .  4436  E  *© 1 

0.39936*01 

0 

0  6  5  50E  o©3 

0 . 48926*01 

0 . 48S9E  *0 1 

0  43016*01 

o 

0  7 8606 ♦OS 

0.27886*01 

0 . 306  2E ♦© t 

0 . 3  2  0  5  E  *0  1 

0 

©  .  9 1 706  *03 

© . 1658E*01 

0  .  1  9  8  9  E  ♦  O  1 

0 . 2  5  4  4  E  *0 1 

© 

0  1 0486*04 

© .  1 709  E*© 1 

0 . 204  7E  *0 1 

0 . 2602E*0 1 

o 

0.11 79E*©4 

©  .  1  74  26*0 1 

©  .  2066  E  *0 1 

0 . 26  42E*0 1 

o 

©  .  13106*04 

© . 1 8006*01 

©  . 2 1 47E*0  1 

© . 2689E*© 1 

0 

©  .  1441  6*04 

©  .  1 8056*0 1 

© . 2 1 57E*0 1 

© .  2  7  O  3  E  *© 1 

0 

0 . 1 5726*04 

0  .  1  8  5  ©  E  *  ©  t 

0 . 2202E*0 1 

0 . 2  7  3  7  E  ♦  0  1 

0 

©  .  1 7036*04 

©  .  1  8526*0  1 

0 . 2 2 © 8  E *0 1 

0 . 2  7  4  4  E  ♦  0  1 

© 

© . 1 8346*04 

©  .  1  8306*0  1 

© . 2 1 906*0 1 

0 . 2  7  3  9  E  *0 1 

0 

©  .  1  1656*04 

© . 1 1016 ♦  © 1 

0 . 21716*01 

0 . 273  2E  *0 1 

0 

© . 20866*04 

0  18516*01 

0 . 221 2E*0 1 

0 . 275  8E  *© 1 

0 

© . 22276*04 

0  .  1  8  7  6  E  ♦  0  1 

0 . 2  2  3  6  E  *  O  1 

0 . 277  3  E ♦© 1 

o 

© . 23586*04 

0  .  18856*01 

©  .  2245E*01 

0 . 21806*01 

0 

©  .  24896*04 

©  .  1  9  2  5  E  ♦  ©  1 

0 . 228  1  E*0 1 

0 . 28016*01 

0 

© . 26206*04 

©  .  1  9096*0 1 

O  .  226  8E*0 1 

0.27956*01 

© 

©  .275  1  E  *-04 

0  .  1  9  406*0 1 

0 . 2  2  9  6  E  *  O  1 

0.28116*01 

0 

©  .  28826*04 

©  .  5  5866*00 

0  .  10496*01 

0.21526*01 

© 

©  .  30 136*04 

0.85676*01 

© . 67006*00 

0 . 204 1 E ♦© 1 

o 

©.31 44E*04 

©  .  2S73E*0  1 

0 . 2709  E  *0 1 

0.30386*01 

© 

0 . 327SE+04 

039246*01 

0.38106*01 

© . 25666*0 1 

o 

0 . 34066*04 

0.38106*01 

0 . 3727E*0 1 

0.35416*01 

© 

©  .  35376*04 

0 . 37251*0  1 

O . 3662E*0 1 

0 . 35201*01 

© 

© . 36686*04 

0 . 38«OE*©1 

0 . 36 1 2E  *0 1 

0.31016*01 

0 

34006*01 

0 . 34006*0 1 

0 . 34006*0 1 

0 . 34006*0 1 

34006*0 1 

0 . 3  4  ©  0  E ♦©  1 

0 . 34006*0 1 

0 . 34006*0 1 

34006*0 1 

0 . 3  400E *© 1 

©  .  34006*0 1 

©  .  34006*0 1 

3  4  OO  E  *0 1 

0 . 34006*0 1 

0 . 34006*01 

© . 34006*0 1 

34076*0 1 

0 . 28  1 56*0 1 

©  .  2257E*0  1 

0.2  1  806*0 1 

346 1 6*0 1 

0.25816*01 

0  1 9076*0 1 

0  16916*01 

33346*0 1 

0.35196*01 

0 . 37976*01 

0.41 776*01 

327  1  E  *  0  1 

040896*01 

0 . 48956*0  1 

0  .  S56  8  E ♦ ©  1 

330  16*01 

0.40656*01 

0.48086*01 

0.5*456*01 

33236*01 

0  40496*01 

0  47506*01 

©  .  S  3  5  8  6  *  0  1 

3338**0 1 

©  .  40206*0 1 

0.46766*01 

0 . 52606*0 1 

334  7E*0 1 

0 . 401  66*01 

0,46576*01 

0 . 52346*0 1 

3  3  S  5  E  *0 1 

0  39936*01 

0.46056*01 

0.51656*01 

33596*01 

O  .  39896*01 

0 . 45966*01 

0.51566*01 

336  4  E*0 1 

0 . 40006*0 1 

0.46086*01 

0.51 706*0 1 

33688*01 

0  40  1  16*01 

0.46236*01 

0 . 51886*01 

23726*01 

0 . 39906*0 1 

0 . 4  S79  6*0 1 

0.51266*01 

33746*0 1 

0 . 39776*0 1 

0 . 45556*0  1 

0 . 50966*01 

33766*01 

0.39736*01 

0 . 4S44E*0  1 

0  50806*01 

33766*0 1 

0,39536*01 

©  .  45096'*©  1 

0 . 50356*01 

33766*01 

0 . 39606*0 1 

0 . 45206*0  1 

0 . 50526*0  1 

33776*01 

0.39456*01 

0.44936*01 

0  50176*01 

3  4  2  2  E  ♦  O  1 

0.46266*01 

©  .  56  506*0  1 

0 . 64436*0 1 

35646*0 1 

0.49146*01 

0.59216*01 

0.64  956*0 1 

343  1  E  *0 1 

0 . 379  6  6*0 1 

0.40606*01 

0.41586*01 

3329  E*0 1 

0.31576*01 

0 . 306  16*0  1 

0 . 303  16*0  1 

33536*0 1 

0 . 32076*0 1 

0.31216*01 

0.30936*01 

33666*9 1 

0 . 32456*0 1 

0.31696*01 

0.31446*01 

33786*0 1 

0 . 32766*0 1 

©  .  32096*0  1 

0.3  1  876*0 1 

Figure  B2.  CRVCNL  Example  Output 
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c 

c  CRVCNL 

C 

C 

C  PROGRAM  TO  CALCULATE  BED  TOPOGRAPHY  OF  A  MEANDERING  RIVER. 

C  USES  METHOD  OF  ENGL EUND (  1 974 )  GENERALIZED  TO  ALLOW 
C  ARBITRARY  MEANDER  GEOMETRY. 

C 

C 

IMPLICIT  REAL *8  (A-H,  0-Z) 

DIMENSION  XM( 100) , YM( 100) ,C( 100) ,DC( 100) 

DIMENSION  DEP1( 101 ,21) ,  U( 101 ,21) ,  DU( 101 ,21) 

DIMENSION  V( 101 ,21 ) ,DUI ( 101 , 21 ) ,DEPC( 101,21) ,DEP2( 101,21) 

INPUT  DATA 

RE AD ( 5  ,  1 00 ) D , B , F , PHI , P 

100  FORMAT ( 5D 10 . 4 ) 

TPHI=7 . *DTAN(PHI*3 . 14159/180. ) 

READ( 5 , 101) NSP , NRP , DS 

101  FORMAT(2I3,D10.4) 

DO  1 0  1  =  1, NSP 

READ ( 5 , 1 02 ) XM ( I ) , YM( I ) 

102  FORMAT ( 2D  10 . 4 ) 

10  CONTINUE 

CALCULATE  CURVATURES 

NSP 1 =NSP- 1 
DO  20  1=2, NSP 1 

C( I ) =DAT  AN( ( YM( 1  +  1 )-YM( I ) )/(XM( 1+1 )-XM( I ) ) ) 

*  -DATAN( ( YM( I )-YM( I- 1 ) )/(XM( I )-XM( I- 1 ) ) ) 

C(I)=C(I)/DS 

IF(XM(I)  .GE.XM(I-I)  . AND . XM ( I )  . GE . XM ( I + 1 ) ) C ( I ) =C ( I  -  1 ) 

I F ( XM ( I ) .LE.XM(I-I)  . AND . XM ( I )  . L E . XM( I  +  1 ) ) C ( I ) =C ( I  -  1 ) 

20  CONTINUE 
DC ( 1 ) =0 . 

DC ( NSP ) =0 . 

C( 1 ) =C ( 2 ) 

C ( NSP ) =C ( NSP 1 ) 

DO  30  1=2, NSP 1 

30  DC( I )  =  (C( 1+ 1 )-C( I- 1 )  )/2 ,/DS 

CALL  OUTPUT ( NSP, DS, 1 , 'CURVATURE' ,9,C) 

CALL  OUTPUT(NSP,DS, 1 , 'DERIVATIVE  OF  CURVATURE ', 23 , DC ) 

CALCULATE  FIRST  APPROXIMATION  TO  BED  TOPOGRAPHY 

DR  =  B/(NRP-  1  ) 

R0=-B/2 . 

D03 1  1=1, NSP 
DO  31  U= 1 , NRP 

3  1  DEP1  (  I  .  J)  =  (  1  .+(  (vJ-1  )  *DR  +  RO )  *C  (  I  )  )**TPHI 

CALL  OUTPUT(NSP,DS, NRP,  'FIRST  DEPTH  APPROX  I  MAT  I  ON '  , 25 , DEP 1 ) 

CALCULATE  LONGITUDINAL  VELOCITY  CORRECTION 

DO  32  J=1 , NRP 
R=( J-1 ) *DR+RO 


B3 «  CRVCNL  Program  Listing 


61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 
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79 
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81 

82 

83 

84 

85 

86 

87 
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U(  1  , J)=R*TPHI*C(  1 )/2. 

DU( 1 , J)=0. 

D032  1=2, NSP 

UK1=DS*(R*(-DC( I- 1 )+TPHI *F/D/2 . *C( I  -  1 ) ) -F/D*U( I  -  1 , J ) ) 

DU( 1-1 ,d)=UK1/DS 

UK2  =DS* ( R*  ( ( -DC ( I  -  1 ) -DC ( I ) )/2 . +TPHI*F/D/4 . * ( C ( I  -  1 ) +C ( I ) ) ) -F/D 
**(U( 1-1 , J)+UK1/2. ) ) 

UK3  =  DS * ( R * ( ( -DC ( 1-1 ) -DC ( I ) )/2 . +TPHI*F/D/4 . * ( C ( I  -  1 ) +C ( I ) ) ) -F/D 
*  * ( U ( I  -  1 , J)+UK2/2 . ) ) 

UK4=DS*(R*(-DC( I )+TPHI*F/D/2 . *C ( I ) ) -F/D* ( U ( I -  1 , J)+UK3) ) 

U( I ,J)=U( 1-1 ,J)+(UK1+2. *UK2+2  .  *UK3+UK4)/6 . 

CONTINUE 

CALL  OUTPUT(NSP,DS,NRP,  'LONGITUDINAL  VELOCITY  CORRECT I  ON '  , 32 , U ) 
CALCULATE  LATERAL  VELOCITY 
D033  1=1, NSP 

T  E  =  (  —  1  . +TPHI ) *DC ( I)+TPHI*F/2./D*C( I ) -F/B/D*U( I ,NRP) 

T I = -DC ( I ) +TPHI *F/2 . /D*C ( I ) - F/B/D*U ( I ,NRP) 

DO  33  d=1 , NRP 
R=( J- 1 ) *DR+RO 

V( I ,U)=TE/2./DEP1 ( I , J)/(  1  ,+C( I )*R) *(RO*RO-R*R) 

DU  I ( I , J)=TI*(R*R-R0*R0)/2 . 

CONTINUE 

CALL  OUTPUT(NSP,DS , NRP,  'LATERAL  INTEGRAL  OF  DU/DS '  , 25 . DU  I ) 

CALL  OUTPUT(NSP,DS, NRP, 'LATERAL  VE LOC I T Y ' , 1 6 , V ) 


Continued 
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90 

C 

CALCULATE  SECOND  APPROXIMATION  TO  DEPTH 

91 

C 

92 

DO  34  1=1, NSP 

93 

T0  =  0 . 

94 

DEPC (1,1) =0 . 

95 

DO  35  J=2 , NRP 

96 

R=( d-1 ) *DR+RO 

97 

T 1  =  (  1  . /DEP 1 ( I,d)*(-1)*(V(I , J)+P/( 1  . +U( I , J) )**P/( 1 . +R*C( I ) ) 

98  ■ 

* *DU I ( I , J) )  ) 

99 

DEPC  (  I\d)  =  (T1+T0)/2  .  *DR*TPHI/7  .  /D  +  DEPC(  I  ,  J-  1  ) 

100 

35 

TO  =  T  1 

101 

NRP 1 =NRP- 1 

102 

DEPK  =  0. 5* (DEP 1 ( I ,  1  )  *DEPC( I  ,  1  )+DEP1 ( I , NRP) *DEPC( I , NRP) ) 

103 

DO  36  U=2 , NRP 1 

104 

36 

DEPK=DEPK+DEP1( I,J)*DEPC(I,U) 

105 

DEPK=( -  1 ) *DEPK/NRP 1 

106 

DEP1A=0.5*(DEP1(I,  1)+DEP1(I, NRP) ) 

107 

D038  J=2 , NRP 1 

108 

38 

DEP1A=DEP1A+DEP1 ( I . J) 

109 

DEP  1  A  =  DEP 1 A/NRP 1 

1  10 

DEPK=DEPK/DEP1A 

1  1  1 

DEPK2= 1  . -DEP  1A 

1  12 

DO  37  J=1 , NRP 

1  13 

37 

DEP2(I,J)=(DEP1(I ,  J )  *  ( 1 .+DEPC(I , J ) +DEPK ) +DEPK2 ) *0 

1  14 

34 

CONTINUE 

1  15 

C 

1  16 

CALL  OUTPUT(NSP,DS, NRP, 'DEPTH  CORR ECT I ON ' , 1 6 , DEPC ) 

m 

CALL  OUTPUT ( NSP , DS , NRP , 'BED  TOPOGRAPHY 1 4 , DEP2 ) 

1  18 

C 

1  19 

STOP 

120 

END 

121 

C 

122 

C 

123 

SUBROUTINE  OUTPUT ( NSP , DS , NS I , T I TL E , LT , VAR ) 

124 

IMPLICIT  REAL*8  (A-H,  0-Z) 

125 

LOGICAL* 1  TITLE ( LT ) 

126 

DIMENSION  VAR ( 101,21) 

127 

WRITE ( 6 , 100) 

128 

100 

FORMAT ( ' 1 ' , 5X , '  CURVED  CHANNEL  BED  TOPOGRAPHY  OUTPUT'/) 

129 

WRITE(6, 101 ) (TITLE(I) , 1=1 ,LT) 

130 

101 

FORMAT (/5X, 'S(M) ' , 3X.80A1//) 

13  1 

WR I TE ( 6 , 105) 

132 

1.05 

FORMAT ( 10X) 

133 

DO  10  1=1, NSP 

134 

S  =  DS*( I-  1  ) 

135 

WR I TE ( 6 , 1 02 ) S , ( VAR ( I , J) , J= 1 . NS I ) 

136 

102 

FORMAT) /I  2D  1  1  . 4/ 1 1 X ,  1  1 D  1  1  . 4 ) 

137 

10 

CONTINUE 

138 

RETURN 

139 

END 

End  of  file 
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PROGRAM  PERTRB 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  V0R(51 ,51) ,C(3,51)  ,RHS(51 ,51) ,PSI(51 .51) ,  C 1 ( 3 . 5  1  ) 
DIMENSION  V(51 ,51 ) ,W(51 ,51 ) ,U(51 ,51) , ET A ( 5 1 ,51 ) 

DIMENSION  DWDY ( 5 1 , 5 1 ) ,  DWDZ(51,51),  UAV(51,51),  D2U0DY (51,51) 
DIMENSION  U0( 51 , 51 ) , DUODY (51,51), DUODZ( 51 , 51 ) ,DU0DYZ(5 1,51) 
DIMENSION  C2(3, 51) ,  UF I NAL ( 1 0 1 , 5 1 ) ,  RHS2(51,51) 

N=5 1 
M  =  5  1 
M0  =  5 


N0  =  5 

MAX  XT  =  200 

OPEN ( UNI T  =  4 ,  FILE3 'DEVDAT ' ,  ST  ATUS= ' OLD ' ) 

OPEN( UNIT=7 ,  FILE3' OUTPUT ' ,  STATUS3 ' NEW ' ) 

CALL  INITA(V0R,M,N,O.O) 

CALL  INITA(PSI .M.N.O.O) 

CALL  INITA(ETA, M.N.O.O) 

READ ( 4 , * ) EDDY ,  ELAT,  CF , BETA , SOR , UOBC , ALPHA 
I F ( EDDY ,LT .0.0)ST0P 
WR I T  E ( 7 ,  '  (  '  '  1  "  )  '  ) 

WR I TE ( 7 , * )  'EDDY  VISCOSITY  COEFFICENT  =  '.EDDY 

WR I TE ( 7 , *  )  'LATERAL  EDDY  VISCOSITY  MULTIPLIER  3  '.ELAT 

WR  I  TE ( 7 , * )  'FRICTION  COEFFICENT  =  '  ,  CF 
WR I TE ( 7  ,  * )  'ASPECT  RATIO  3  ',  BETA 

WR I TE ( 7 , * )  'SUCCESIVE  OVER-RELAXATION  FACTOR  3  ',  SOR 
WR I T  E  (  7  , *  )  'SIDE  BOUNDARY  CONSTANT  =  '  , UOBC 
WR I T E ( 7 , * )  'WIDTH  TO  RADIUS  RATIO  3  '.ALPHA 
CALL  INITA(UO, M.N.O.O) 

SOR  =  SOR/( 1  .0+1  . O/BETA  *EL AT ) 

UBB  3 1 .0-1 .0/3 .O/CF/EDDY 
UBS3 1 .0+1 .0/6 . O/CF/EDDY 
UORHS  3 -  1  . O/EDDY/CF 
CALL  INITA(RHS,M,N, UORHS) 

DO  33  1=1 ,M 
C(  1  ,  I)  =  1 
Cl  (  1  ,  I  )  =  1 

Cl  (2,  I )  =  1  .  0/  EDDY/ CF/UBB  *B  ETA/ EL AT  *UOBC 
Cl (3, I) =0.0 
C2 (  1  , I ) 3 1 

C2(2, I)=1 . O/EDDY/CF /UBB *BET A/ ELAT 
C2 ( 3 , I ) =0 . 0 

C ( 2  ,  I ) 3 -  1  . O/EDDY/ CF/UBB 
C ( 3 , I ) =0 . 0 

CALL  I TER A ( UO , M , N ,  1  ,  1 ,2,2,0.00001 ,C1 , C , SOR , BETA , RHS , ELAT ) 
CALL  D I NT EG( UO , M , N , RHS , 2 , 1 .0) 

CALL  DINTEG(RHS,M,N,UAV, 1 , 1 .0) 

DO  15  1=1,  M 
DO  15  d=  1  ,  N 

UAV ( I , J)=UO( I , J)/UAV(M,N) 

CALL  OUTA(UAV, M,N, MO, NO,  ' INITIAL  VELOCITY  D I  STR IBUT ION '  , 29 ) 
CALL  DERIVA(UAV,M,N, DUODZ , 1 , 2, 1 .0) 

CALL  DERIVA(UAV,M,N, DUODY , 1 , 1 , 1 .0) 

CALL  DERIVA(UAV,M,N, D2U0DY , 2 , 1,1.0) 

CALL  DER I VA ( DUODY , M , N , DUOD Y Z ,  1 ,  1 .  1  .0) 

DO  10  J=1 ,N 
DO  10  I3  1  ,M 

RHS  (  I , J)3CF/EDDY*2. *UAV( I . J)*DUODZ( I . J) 

CALL  ITERA(VOR,N,N, 1,0, 0,0.0. 00001 . C 1 . C , SOR , BETA , RHS . ELAT ) 
DO  22  J=1,N 
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6  1 

00  22  I = 1 , M 

62 

22 

RHS2( I , J) =-V0R( I , J)/B ETA/BETA 

63 

CALL  ITERA(PSI ,N,N,  1 ,0,0,0,0.00001 ,C1 , C , 50R , BETA , RHS2  ,  1.0) 

64 

VC=BETA 

65 

CALL  DERIVA(PSI , N , N , V ,  1  ,2, VC) 

66 

CALL  DERIVA(PSI ,N,N,W,  1 ,  1 , -  1  .0) 

67 

I F ( ABS ( V ( 1 , 1 ) -VO) . LT. 0.001 )GOTO  31 

68 

VO=V( 1,1) 

69 

DO  30  1=1 ,M-1 

70 

30 

VOR( I ,  1  )  =-BETA/EDDY/CF*V( I ,  1  )/UBB 

71 

DO  32  U  =  2,N 

72 

32 

VOR(M,  J)  =-W(M,  J)/EDDY/CF/UBB*-BETA 

73 

GOTO  21 

74 

3  1 

CONTINUE 

75 

CALL  OUTA(V,M,N,MO,NO, 'LATERAL  VELOCITY  D I STR I BUT I ON '  , 29 ) 

76 

CALL  OUTA(W,M,N, MO, NO, 'VERTICAL  VELOCITY  D I STR I BUT I ON ', 30 ) 

77 

CALL  D I NTEG ( W,M,N,RHS2,2, 1 .0) 

78 

CALL  OUT A ( RHS2 , M , N , MO , NO , 'DEPTH  AVG  VERTICAL  VELOC I TY ' , 27 ) 

79 

CALL  DERIVA(PSI ,M,N,DWDY,2, 1 , -1 .0) 

80 

CALL  DERIVA(W,M,N,DWDZ. 1 ,2. 1 .0) 

8  1 

DO  40  1=1 ,M 

82 

DO  40  J=1 ,N 

83 

RHS( I  , U ) = - ( V ( I , J) *BETA*CF/EDDY+ 1 ) *D2U0DY( I , J ) /BETA/BET  A 

84 

RHS ( I , U ) =RHS ( I , J)-CF/EDDY*(DWDY( I ,  U)*DUODZ( I , J) 

85 

*  -DWDZ ( I , J)*DUODY( I , J ) +W( I ,U)*DUODYZ( I ,  U  )  +  1 .O/CF/CF) 

86 

40 

CONTINUE 

87 

I  T  =  0 

88 

US0=-9999 . 9 

89 

42 

I T  =  I T+ 1 

90 

C 

IF( IT.GT ,MAXIT)GOTO  43 

91 

CALL  I TERA ( ET  A , M , N ,  1 ,  1 ,0,2,0.0001 ,C1 , C , SOR , BETA , RHS , EL AT ) 

92 

CALL  DINTEG(ETA,M,N,U, 1 ,-1 .0) 

93 

IF(ABS(U(M,N)-USO) . LT . 0 . 002 ) GOTO  50 

94 

USO=U ( M , N ) 

95 

WRITE ( *  ,  *  )  '  USO  =  '  ,  USO 

96 

DO  51  J=  1  ,N 

97 

51 

ETA(M, J) =U(M, J)/EDDY/CF/UBB*BETA/ELAT*UOBC 

98 

GOTO  42 

99 

50 

CONTINUE 

100 

CALL  OUTA(U,M,N, MO, NO, 'LONGITUDINAL  VELOCITY  D I STR I BUT I ON ', 34 ) 

101 

CALL  D I NTEG ( U,M,N,RHS,2, 1 .0) 

102 

CALL  OUTA(RHS,M,N, MO, NO, 'DEPTH-AVERAGED  LONG.  VELOC I T Y ' , 29 ) 

103 

DO  6  1  1  =  1  ,N 

104 

UF INAL ( M , I ) =UA V (  1  , I  )+ ALPHA *U(  1,1) 

105 

DO  61  J=1 , M- 1 

106 

UP=M+J 

107 

UM=M- J 

108 

UF INAL (UP, I)=UAV( J+1 , I  ) + AL PHA *U ( U+ 1  ,  I) 

109 

UF I NAL ( UM , I )=UAV(J+1 , I ) -ALPHA*U(  J+  1,1) 

1  10 

6  1 

CONTINUE 

1  1  1 

MF=2*M- 1 

1  12 

CALL  OUTA ( UF I NAL , MF , N , MO , NO ,  'TOTAL  VELOCITY  D I STR I BUT I ON ', 27 ) 

1  13 

GOTO  1 

1  14 

C 

43 

WRITE ( 7 , * )  '  MAXIMUM  NUMBER  OF  SIDE  BOUNDARY  ITERATIONS' 

1  15 

C 

GOTO  1 

1  16 

END 

1  17 

C 

1  18 

C 

1  19 

SUBROUTINE  ITERA(A ,M,N, NB  C 1 , NBC 2 , NBC 3 , NBC 4 ,T0L,C2,C,W,B,RHS,E) 

120 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 
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DIMENSION  A(M,N),C(3,*), RHS(M, N) ,C0( 3 , 200) , C 1(3,200) . C2( 3 , 200) 
A0= -999999 . 999 
A  1  =  A0 

MAXI T  =  6000 
I T  =  1 

M2  = ( M- 1 )/2 

N2=(N- 1 )/2 

H  = 1 . 0/B/B*E 

DXD Y  = 1 , 0/ (M-1 )/ (N-1 ) 

R=0. 25*W 

CC1=1 . 0-2.0*R*H-2.0*R 

CC2=R*H 

CC3  =  R 

CC4  =  - R  *DXD Y 
DO  20  I = 1 , M 
C0( 1 . I)=1 .0 
C0( 2 , I)=0.0 

20  C0( 3 , I ) =0 . 0 
DO  21  1=1 ,N 
Cl  (  1  ,  I  )  =  1 .0 
C 1 ( 2 , I)=0.0 

21  C 1 ( 3 , I ) =0 . 0 

10  CALL  AASOR(A,M,N,W,B,RHS,E) 

I F ( NBC  1  . EO . 2 ) CALL  BCSOR(A,M,N,  1 , W , B , RHS , C , E ) 

I F ( NBC  1  . EO.  1  ) CALL  BCSOR(A,M,N,  1 , W . B . RHS , C 1 ,E) 

IF(NBC2 . EO . 2) CALL  BCSOR ( A , M , N , 2 , W , B , RHS , C , E ) 

I F ( NBC2 . EO .  1)CALL  BCSOR ( A , M , N , 2 , W , B , RHS , CO . E ) 

IF(NBC3 . EO . 2) CALL  BCSOR ( A , M , N , 3 , W . B , RHS , C2 , E ) 

I F ( NBC3 . EO .  1)CALL  BCSOR ( A , M , N , 3 . W , B , RHS , C 1 , E) 

I F ( NBC4 . EO . 2) CALL  BCSOR ( A , M , N , 4 , W, B , RHS , C , E ) 

I F ( NBC4 . EO .  1  ) CALL  BCSOR ( A , M , N , 4 , W , B , RHS , CO , E ) 

I F ( NBC3 . EO. 2 . AND . NBC4 . EO . 2 ) THEN 

AI=A(M-1 ,  1)-C2(2,  1  ) * A ( M ,  1)*2. 0/ ( M- 1 ) 

A J= A ( M , 2 ) +C ( 2 ,M)*A(M, 1  )  *  2  . 0/ ( N- 1 ) 

A ( M , 1 )=CC 1  *A(M,  1 )+CC2*( AI+A(M- 1,1) ) +CC3 * ( A U+ A ( M , 2 ) ) 

*  +CC4  *  RHS ( M ,  1 ) 

ENDIF 
I T= I T+  1 

IF(  IT .GT .MAXIT )GOTO  31 

IF(ABS( (A(M2,N2)-A0)/A0)  . L T . TOL ) GOTO  32 
A0= A ( M2 , N2 ) 

GOTO  10 

30  IF(ABS( ( A ( M- 1 , N- 1 )-A1 ) / A  1 )  . LT. TOL) GO TO  32 
A  1 = A ( M-  1  , N- 1 ) 

GOTO  10 

31  WRITE(*,*)  '  MAXIMUM  NUMBER  OF  ITERATIONS  REACHED' 

32  RETURN 
END 


SUBROUTINE  INI TA ( A , M , N . C ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.  0-Z) 
DIMENSION  A ( M , N ) 

DO  10  1= 1 ,N 
DO  10  J= 1 ,M 
10  A ( J , I  )  =C 
RETURN 
END 
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SUBROUTINE  AASOR(A,M,N,W,B, RHS , E ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.  0-Z) 

DIMENSION  A ( M , N ) , RHS ( M , N ) 

H= 1 . 0/B/B*E 

DXDY  =  1  . 0/ ( M- 1 ) / ( M- 1 ) 

R=0 . 25*W 

C  1  =  1  0_2 .  0  *  R  *  H  -  2 ,0*R 

C2=R*H 

C3  =  R 

C4  =  - R  *DXDY 
DO  10  J=2 , M- 1 
DO  10  K=2,N-1 

10  A(U,K)=C1*A(J,K)+C2*(A(J- 1 , K )  +  A  (  J+ 1  , K ) ) +  C3* ( A ( J , K- 1 ) +  A ( d , K+ 1 ) ) 
*+C4*  RHS ( J , K ) 

RETURN 

END 

C 

C 

SUBROUTINE  BCSOR ( A , M , N , NB , W , B , RHS , C , E) 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  A ( M , N ) ,RHS(M,N) ,C(3, *  ) 

H= 1 . 0/B/B*E 

DXD Y  =  1  . 0/ ( M- 1 ) / ( N- 1 ) 

R=Q , 25  *W 

C 1 = 1 .0-2.0*R*H-2.0*R 

C2=R*H 

C3  =  R 

C4  =  -R *DXDY 
I F ( NB . EQ. 1 ) THEN 
DO  10  1=2, N-1 

AI=A(2, I )-2 .0*(-C(2, I )*A( 1 , I )+C(3, I ) )/C( 1 , I )/(M-1 ) 

10  A(  1  , I ) =C 1 *A(  1 , I  )+C2*( AI  +  A( 2 ,  I  )  )+C3*( A(  1  ,  I  -  1  )  + A (  1  , I +  1  ) ) 

*  +C4*RHS (1,1) 

ELSEIF(NB . EO. 2)THEN 

AI=A( 1 ,N-1 )+2 .0*(-C(2 ,  1  )*A(  1  ,N)+C(3,  1  )  )/C( 1 ,  1  ) / (N- 1 ) 

A ( 1 , N) =C 1 *A ( 1 , N)+C3*( AI+A( 1 , N- 1 ) ) +C2* ( A ( 2 , N ) + A ( 2 , N) ) 

*  +C4*RHS(1.N) 

DO  20  1=2, M-1 

A  I =A( I , N- 1 )+2 .  0*( -C( 2 , I ) *A( I , N)+C( 3 , I ) )/C( 1 , I )/(N- 1 ) 

20  A(I,N)=C1*A(I , N ) +C3  * (  A  I  +  A  ( I , N- 1 ) )+C2*( A(  I  -  1  ,  N  )  +  A  (  I  +  1  .  N )  ) 

*  +C4*RHS ( I , N ) 

ELSEIF(NB.E0.3)THEN 

DO  30  1=2, N-1 

AI=A(M-1 , I )+2.0*(-C(2, I )*A(M, I )+C(3, I ) )/C( 1 , I )/(M-1 ) 

30  A(M,I)=C1*A(M,I ) +C2  * ( A I + A ( M- 1,1) ) +C3* ( A ( M , I  -  1 ) + A ( M , I  +  1 ) ) 

*  +C4*RHS(M,I) 

AI=A(M- 1 ,N)+2 ,0*(-C(2,N)*A(M,N)+C(3,N) )/C( 1 ,N)/(M- 1 ) 

A ( M , N ) =C 1 *A(M,N)+C2*(AI+A(M- 1 , N ) ) +C3* ( A ( M , N- 1 )+A(M, N- 1) ) 

*  +C4  *RHS ( M , N ) 

ELSE  I F ( NB . EO . 4 ) THEN 

AI=A( 1 ,2) -2 .0*(-C(2, 1 )*A( 1 , 1 )+C(3, 1 ) )/C(  1  ,  1 )/(M-1 ) 

A( 1  ,  1 )=C1*A(  1  ,  1 )+C2*( A(2  ,  1  )+A(2,  1  ) )+C3*( AI  +  A(  1,2)) 

*  +C4*RHS( 1,1) 

DO  40  1=2, M-1 

AI=A(I,2)-2.0*(-C(2,I)*A(I,1)+C(3,I))/C(1,I)/(M-1) 

40  A ( I ,  1)=C1*A(I,  1  )+C2*( A( 1-1  .  1  )+A( 1+1  ,  1  ) ) +C3 * ( A  I +A ( I , 2)  ) 

*  +C4*RHS (1,1) 

ELSE 

WRITE ( * , 200 )NB 

200  F0RMAT(5X , 'THE  VALUE  OF  NB  IS  INCORRECT  ',13/) 
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24  1 

END I  F 

242 

RETURN 

243 

END 

244 

C 

245 

C 

246 

SUBROUTINE  OUT A ( A , M , N , NX , NY , T I TL E 

247 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0 

248 

DIMENSION  A ( M , N ) , TOP (21) 

249 

CHARACTER  TITLE(NCT) 

250 

WR I T  E ( 7 , 200 ) ( T I TL  E ( I ) ,1  =  1 ,NCT) 

251 

200 

FORMAT ( '0' , 80A 1 // ) 

252 

DX= 1 . 0/(M- 1 ) 

253 

D Y  = 1 . 0/ ( N- 1 ) 

254 

N0= ( N- 2 ) /NY 

255 

M0= ( M-2 ) /NX 

256 

DO  10  1=1, MO 

257 

10 

TOP ( I ) = I *NX*DX 

258 

TOP ( M0+  1  )  =  1  .0 

259 

WRITE(7, 201 ) ( TOP ( I ) , 1= 1 , M0+ 1 ) 

260 

201 

FORMAT ( ' 0 ' , 1 OX , '  0.0  '.21F10 

26  1 

DO  20  d=N, NY+ 1 , -NY 

262 

K  =  0 

263 

DO  30  1=1, M-NX , NX 

264 

K  =  K+  1 

265 

30 

TOP ( K ) =A ( I , U) 

266 

Y  = ( U  —  1 )  *DY 

267 

20 

WR I TE ( 7 , 202 ) Y , (TOP(K) ,K=1 , M0+ 1 ) .A 

268 

202 

FORMAT ( IX , 22F 10. 4) 

269 

K  =  0 

270 

DO  40  1  =  1  , M-NX, NX 

27  1 

K  =  K+  1 

272 

40 

TOP(K)=A( 1,1) 

273 

Y  =0 . 0 

274 

WR ITE ( 7 , 202 ) Y , (TOP(K) , K=1 ,M0+1) , A i 

275 

RETURN 

276 

END 

277 

C 

278 

C 

279 

SUBROUTINE  ADUMP ( A , M , N ) 

280 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0- 

28  1 

D I  MENS  I ON  A ( M , N )  ,  0UTV(13) 

282 

MR  =  M 

283 

MIN=  1 

284 

MAX=  1  3 

285 

WR I T E ( 7 , * )  '  MATRIX  DUMP' 

286 

30 

MO  =  MR 

287 

I F ( MR . GT . MAX ) MO=MAX 

288 

DO  10  1=1 ,N 

289 

DO  11  J= 1 , MO 

290 

1  1 

OUTV(J)=A( MI N+ J- 1 . I) 

291 

WR I TE ( 7 , 200) (OUTV(J) , U= 1 , MO ) 

292 

200 

FORMAT ( IX, 13F 10. 4) 

293 

10 

CONTINUE 

294 

I F ( MO . EO . MR ) GOTO  20 

295 

MR=MR-MAX 

296 

MIN=MAX+MIN 

297 

WR ITE ( 7 , *  )  '  MATRIX  CONTINUATION' 

298 

GOTO  30 

299 

20 

RETURN 

300 

END 
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SUBROUTINE  DERIVA(A,M,N,D, ND , NYZ , C ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  A ( M , N ) ,  D ( M , N ) 

DX  = 1 . 0/ ( M- 1 ) 

DY  = 1 . 0/ ( N- 1 ) 

I F ( ND . EQ . 1)  THEN 
I F ( NY  Z . EQ.  1 )  THEN 
DO  10  I  =  1  »  N 

D(1,I)=(A(2,I)-A(1,I))/DX 
DO  11  J=2 , M- 1 

D(J.I)=(A(J+1,I)-A(J-1,I ) )/2 ./DX 
D(M, I ) = ( A ( M , I ) -A ( M- 1 ,1) ) /DX 
ELSEIF(NYZ. EQ . 2 ) THEN 
DO  20  1  =  1  ,M 

D (  I  ,  1)=(A(I,2)-A(I,  1))/DY 
DO  21  J=2 , N- 1 

D( I , J)  =  ( A( I , J+1  ) -A( I ,J- 1 ) )/2 . /DY 
D(I,N)=(A(I,N)-A(I,N-1 ) )/DY 
ELSE 

WR I  T  E '  THE  VALUE  OF  ND  IS  OK  BUT  NOT  NX Y '  ' ) ' ) 

END  I F 

ELSE  I F ( ND . EO . 2 )  THEN 
I F ( NYZ . EO . 1 )  THEN 
DO  30  1=1 .N 

D(1,I)=(2.*A(2,I)-2.*A( 1,1 ) )/DX/DX 
DO  31  J=2,M-1 

D(J, I)=(A(J+1  , I  )  -  2  . *  A  (  J  , I ) + A ( J- 1 , I ) ) /DX/DX 

D(M, I )=(2. *A(M, I )-5 . *A(M-1 , I )+4 . *A(M-2, I )-A(M-3 . I ) ) /DX/DX 
ELSEIF(NYZ . EO . 2 ) THEN 
DO  40  1=1 ,M 

D ( I ,  1)  =  (2.*A(I,2)-2.*A(I,  1  )  )  /DY/DY 
DO  41  J=2 , N- 1 

D( I . J)=(A( I , d+1)-2. *A( I , J)+A( I , J- 1 ) ) /DY/DY 

D( I ,N)=(2. *  A  ( I ,N)-5. *A( I ,N-1 )+4. *A( I , N-2 ) -A ( I ,N-3) ) /DY/DY 
ELSE 

WR I TE '  THE  VALUE  OF  ND  IS  OK  BUT  NOT  NXY '  ' ) ' ) 

END  IF 
ELSE 

WR I TE '  THE  VALUE  OF  ND  IS  I NCORRECT '  ' ) ' ) 

END  IF 

DO  50  1= 1 ,N 
DO  50  J= 1 , M 
D( J, I )=C*D( J, I ) 

RETURN 

END 


SUBROUTINE  DINTEG(A,M,N,DIN, NY Z . C ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  A ( M , N ) ,  D I N ( M , N ) 

DX  = 1  . 0/ ( M- 1 ) 

DY= 1 . 0/ ( N- 1 ) 

CALL  INITA(DIN,M,N,0.0) 

I F ( NYZ . EO. 1 )  THEN 
DO  10  U= 1 . N 
DO  10  1=2, M 

DIN( I , J)=DIN( I- 1 , d)+0 . 5*( A(I- 1 ,J)+A( I . J) )*DX*C 
ELSEIF(NYZ . EO. 2)  THEN 
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36  1 

DO  20  I=1,M 

362 

DO  20  0  =  2, N 

363 

20  DIN( I , J)=DIN( I , 0- 1 )+0. 5*DY*( A( I  ,  0- 1 ) + A ( I , 0) ) *C 

364 

ELSE 

365 

WRITE(*,*)  '  THE  VALUE  OF  NX Y  IS  INVALID  ' 

366 

END  I  F 

367 

RETURN 

368 

END 

of  file 
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2 

3 

4 

5 

6 

7 

8 


PROGRAM  DEVFIO 

IMPLICIT  DOUBLE  PRECISION  (A-H,  O-Z) 

DIMENSION  VOR(21 ,21) ,C(3,21) ,RHS(21 ,21) ,PSI(21 ,21) ,  0(3,21) 

V ( 2 1 , 2 1 ) ,W(21 ,21) ,U(21 ,21) , ET  A ( 2  1,21), ETAT ( 2  1,21) 
DWDY (21,21),  DWDZ (21,21),  UAV(21,21),  D2U0DY (21,21) 
U0( 2 1 , 21 ) , DUOD Y ( 2 1,21) ,DUODZ( 21 , 21 ) . DUODYZ( 21 ,21) 


DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 


C2 ( 3 , 2  1  )  ,  UFINAL(  101 ,21 ) ,  RHS2(21,21),  CDO(3,21) 


9 

N  =  2  1 

10 

M  =  2  1 

1  1 

MO  =  2 

12 

N0  =  2 

13 

MAXI T  =  200 

14 

0PEN(UNIT  =  4 ,  F I L  E  = ' DEVDAT '  ,  ST  ATUS= ' OLD ' ) 

15 

OPEN( UNIT=7 ,  FILE* 'OUTPUT ' ,  ST ATUS = ' NE W ' ) 

16 

CALL  INITA(V0R,M,N,0.0) 

17 

CALL  INITA(PSI .M.N.O.O) 

18 

CALL  INITA(ETA,M,N,0.0) 

19 

1 

RE AD ( 4 , * ) EDDY ,  ELAT,  CF , BETA , SOR . UOBC , ALPHA 

20 

I F ( EDDY . LT . 0 . 0) STOP 

21 

WRITE) 7,  '  (  '  '  1  '  '  )  '  ) 

22 

WR I T  E  (  7  ,  *  )  'EDDY  VISCOSITY  COEFFICENT  =  '.EDDY 

23 

WR I TE ( 7 , * )  'LATERAL  EDDY  VISCOSITY  MULTIPLIER  =  '.ELAT 

24 

WRITE ( 7 , * )  'FRICTION  COEFFICENT  =  ',  CF 

25 

WR I TE ( 7  ,  *  )  'ASPECT  RATIO  =  ',  BETA 

26 

WR I  TE ( 7 , * )  'SUCCESIVE  OVER-RELAXATION  FACTOR  =  ',  SOR 

27 

WRITE)?,*)  'SIDE  BOUNDARY  CONSTANT  =  '.UOBC 

28 

WR I T  E ( 7 , *  )  'WIDTH  TO  RADIUS  RATIO  =  '.ALPHA 

29 

CALL  INITA(UO, M.N.O.O) 

30 

SOR  =  SOR/(  1  .0+1 . O/BET  A  *  EL AT ) 

31 

UBB  = 1 .0-1 .0/3 .O/CF/EDDY 

32 

UBS= 1 .0+1 .0/6. O/CF/EDDY 

33 

UORHS=- 1 . O/EDDY/CF 

34 

CALL  INITA(RHS,M,N, UORHS ) 

35 

DO  33  1  =  1  ,M 

36 

C(  1  , I)  =  1  .0 

37 

C1( 1 , I)  =  1  .0 

38 

Cl (2, I  )  =  1  . 0/ EDDY/ CF /UBB* BET  A/ EL AT  *UOBC 

39 

Cl (3,1) =0.0 

40 

C2(  1  , I  )  =  1  .0 

41 

C2( 2 , I ) = 1 . 0/EDDY/CF/UBB*BETA/ELAT 

42 

C2 ( 3 , I)=0.0 

43 

CDO( 1 , I )  =  1 . 0 

44 

CDO ( 2 , I )=0.0 

45 

CDO( 3 , I )=0.0 

46 

C ( 2 , I ) =- 1 . O/EDDY/CF/UBB 

47 

33 

C ( 3 , I ) =0 . 0 

48 

CALL  ITERA(UO, M , N , 1 , 1 ,2.2,0.00001 ,C1 , C , SOR , BETA , RHS . ELAT ) 

49 

CALL  D I NTEG ( UO ,M,N,RHS,2, 1 .0) 

50 

CALL  DINTEG(RHS,M,N,UAV, 1 , 1 .0) 

5  1 

DO  15  1=1,  M 

52 

DO  15  J=  1  ,  N 

53 

15 

UAV ( I , J)=UO( I . J)/UAV(M,N) 

54 

CALL  OUTA(UAV, M,N, MO. NO,  ' INITIAL  VELOCITY  D I STR I BUT I ON '  , 29 ) 

55 

CALL  DERIVA(UAV,M,N, DUODZ , 1 ,2, 1 .0) 

56 

CALL  DER I VA( UAV, M,N, DUOD Y, 1 , 1 , 1 .0) 

57 

CALL  DERIVA(UAV,M,N, D2U0D Y , 2 , 1,1.0) 

58 

CALL  DER I VA ( DUOD Y ,M,N, DUOD YZ, 1 , 1 , 1 .0) 

59 

DO  10  d=1,N 

60 

DO  10  1= 1 ,M 
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6  1 

10 

RHS( I . J) =CF/EDDY*2 . *UAV( I ,  J) *DU0DZ( I , J) 

62 

2  1 

CALL  I TE  R A ( VOR , N , N ,  1 ,0,0,0,0.00001 , Cl  , C , SOR , BETA , RHS , EL AT ) 

63 

DO  22  J=1,N 

64 

DO  22  1  =  1  ,M 

65 

22 

RHS2 ( I , J ) = -VOR ( I , J) /BETA/BETA 

66 

CALL  ITERA(PSI ,  N  ,  N  , 1 ,0,0,0,0.00001 ,C1 , C , SOR , BETA , RHS2 , 1.0) 

67 

VC  =  BET  A 

68 

CALL  DERIVA(PSI ,  N  ,  N  ,  V  , 1 ,2, VC) 

69 

CALL  DERIVA(PSI ,N,N,W, 1 , 1 , -1 .0) 

70 

IF(ABS(V(  1  ,  1)-V0)  . LT. 0.001 )GOTO  31 

7  1 

VO=V( 1,1) 

72 

DO  30  1=1 , M- 1 

73 

30 

VOR ( I  ,  1 )=-BETA/EDDY/CF*V( I ,  1 )/UBB 

74 

DO  32  J  =  2  ,  N 

75 

32 

VOR(M,J)=-W(M,J)/EDDY/CF/UBB*BETA 

76 

GOTO  21 

77 

31 

CONTINUE 

78 

CALL  OUTA(V,M,N, MO, NO,  'LATERAL  VELOCITY  D I STR I  BUT  I  ON '  , 29 ) 

79 

CALL  OUTA(W,M,N, MO, NO,  'VERTICAL  VELOCITY  D I STR I BUT I ON ', 30 ) 

80 

CALL  D I NTEG ( W,M,N,RHS2,2, 1 ,0) 

8  1 

CALL  OUT A ( RHS2 , M , N , MO , NO , 'DEPTH  AVG  VERTICAL  VELOC I TY ' , 27 ) 

82 

CALL  DERIVA(PSI ,M,N,DWDY,2. 1 ,-1 .0) 

83 

CALL  DERIVA(W,M,N,DWDZ, 1 ,2, 1 .0) 

84 

DX=0. 25*CF/BETA/EDDY/ (M-1 )/ (M-1 ) 

85 

XL  =0 .  1 *CF/BETA 

86 

XE  =  32  ,313 

87 

X=-5 .0 

88 

NSTEP=0 

89 

35 

X=X+DX 

90 

DO  40  1=1 ,M 

91 

DO  40  J=1 ,N 

92 

RHS ( I,J)=-(V(I,U) ) *D2U0DY  (  I  ,  J) 

93 

RHS ( I , J)=RHS( I , J)-BETA*(DWDY( I , J ) *DU0DZ ( I , J) 

94 

*  -DWDZ ( I , U ) *DUOD Y ( I , J)+W( I , J ) *DU0DYZ ( I , d) ) 

95 

RHS ( I , J)=RHS( I , J)*SECF(X,XE,XL)-SUPF(X,XE) 

96 

*-RFUNC(X,XE)*BETA/CF/CF 

97 

40 

CONTINUE 

98 

CALL  AAEXP(ETA, ETAT.M.N, 1 . 0 , BETA , RHS , ELAT , DX ) 

99 

CALL  BCEXP( ETA, ETAT.M.N, 1 , 1 . 0 , BETA , RHS , CDO , ELAT , DX ) 

100 

CALL  BCEXP( ETA , ETAT , M , N , 2 , 1 .0, BETA , RHS , CDO, ELAT, DX) 

101 

CALL  BCEXP(ETA, ETAT,M,N,4, 1 . 0 , BETA , RHS , C , ELAT , DX ) 

102 

NSTEP=NSTEP+ 1 

103 

I F ( NSTEP . EO . 10) GOTO  50 

104 

GOTO  35 

105 

50 

CONTINUE 

106 

CALL  DINTEG(ETA,M.N.U, 1,-1 .0) 

107 

CALL  DINTEG( U,M,N,RHS,2, 1 .0) 

108 

WR I T  E ( 7 ,  121 )X,RHS( 17,N) ,RHS( 19.N) ,RHS(21 ,N) 

109 

121 

FORMAT) 5X , 4F 15 . 5) 

1  10 

NSTEP=0 

1  1  1 

I F ( X . GT .40.0) GOTO 7 1 

1  12 

GOTO  35 

1  13 

71 

GOTO  1 

1  14 

END 

1  15 

C 

1  16 

C 

1  17 

SUBROUTINE  ITERA ( A , M , N , NBC  1 , NBC2 , NBC3 , NBC4 , TOL , C2 , C , W , B , RHS , E) 

1  18 

IMPLICIT  DOUBLE  PRECISION  (A-H.  0-Z) 

1  19 

DIMENSION  A(M.N) , C(3, *) .RHS(M.N) , C0( 3 , 200) , C 1 ( 3 , 200) . C2 ( 3 , 200) 

120 

A0= -999999 . 999 
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121 

A  1  =  AO 

122 

MAX  I T  =  6000 

123 

I T  =  1 

124 

M2= ( M- 1 ) / 2 

125 

N2=(N- 1 )/2 

126 

H= 1 . 0/B/B*E 

127 

DXD  Y  = 1  . 0/ ( M- 1 ) / ( N— 1 ) 

128 

R=0 . 25*W 

129 

CC 1 = 1 . 0-2 . 0*R*H-2 . 0*R 

130 

CC2=R*H 

131 

CC3  =  R 

132 

CC4=-R*DXDY 

133 

DO  20  1  =  1  ,M 

134 

C0( 1 , I )=1 .0 

135 

C0(2, I) =0.0 

136 

20 

C0( 3 , I)=0.0 

137 

DO  21  1  =  1  ,N 

138 

C 1 (  1  ,  I)  =  1  .0 

139 

C 1 ( 2 , I ) =0 . 0 

140 

21 

C 1 ( 3 , I)=0.0 

14  1 

10 

CALL  AASOR(A,M,N,W,B,RHS,E) 

142 

I F ( NBC  1  . EO . 2 ) CALL  BCSOR(A,M,N,  1 , W , B , RHS , C , E ) 

143 

I F ( NBC  1  . EO  -  1  ) CALL  BCSOR ( A , M, N, 1 , W. B . RHS , C 1 , E) 

144 

I F ( NBC2 . EO . 2 ) CALL  BCSOR ( A . M , N . 2 , W . B , RHS , C , E ) 

145 

I F ( NBC2 . EO .  1  )  CALL  BCSOR ( A , M , N , 2 , W , B , RHS , CO , E ) 

146 

IF(NBC3 . EO . 2)CALL  BCSOR ( A . M , N , 3 , W . B , RHS , C2 , E ) 

147 

IF (NBC3 . EO .  1  )  CALL  BCSOR ( A , M , N , 3 , W,B , RHS , Cl , E ) 

148  • 

I F ( NBC 4 . EO . 2) CALL  BCSOR ( A , M . N , 4 , W, B , RHS , C, E) 

149 

IF (NBC4 . EO .  1  )  CAL  L  BCSOR ( A , M . N , 4 , W, B , RHS .CO, E ) 

150 

I F ( NBC3 . EO . 2 . AND .NBC4 . EO . 2 ) THEN 

151 

A  I =A(M- 1 ,  1  )  -C2( 2 ,  1  ) *A(M,  1 )*2 ,0/(M- 1 ) 

152 

AU= A ( M , 2)+C(2,M)*A(M, 1 ) *2 . 0/(N- 1 ) 

153 

A ( M , 1 ) =CC1*A(M, 1 )+CC2*( AI+A(M-1 , 1 ) ) +CC3* ( AU+A ( M , 2 ) ) 

154 

*  +  CC4  *  RHS ( M ,  1 ) 

155 

END  I  F 

156 

I T= IT+ 1 

157 

I  F ( I T . GT . MAX  I T ) GOTO  31 

158 

IF(ABS( (A(M2,N2)-AO)/AO) ,LT.TOL)GOTO  32 

159 

A0= A ( M2 , N2 ) 

160 

GOTO  10 

161 

C 

30 

I F ( ABS ( (A(M-1,N-1)-A1)/A1).LT. TOL ) GOTO  32 

162 

C 

A  1 = A ( M- 1 ,N-1) 

163 

C 

GOTO  10 

164 

31 

WR I TE ( * , *  )  '  MAXIMUM  NUMBER  OF  ITERATIONS  REACHED' 

165 

32 

RETURN 

166 

END 

167 

C 

168 

C 

169 

SUBROUTINE  I NITA ( A , M , N , C ) 

170 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

171 

DIMENSION  A ( M , N ) 

172 

DO  10  1=1 ,N 

173 

DO  10  U=  1  ,  M 

174 

10 

A ( J , I ) =C 

175 

RETURN 

176 

END 

177 

C 

178 

C 

179 

SUBROUTINE  A A SOR ( A , M , N , W . B . RHS . E ) 

180 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 
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18  1 

182 

183 

184 

185 

186 

187 

188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209 

210 

2  1  1 

2  12 

213 

214 

215 

216 

217 

218 

2  19 

220 

221 

222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 

240 


254 


DIMENSION  A(M,N) ,RHS(M,N) 

H= 1 . 0/B/B*E 

DXDY=  1  0/(M- 1 )/(M- 1 ) 

R=0. 25  *W 

C 1 = 1 .0-2.0*R*H-2.0*R 

C2=R*H 

C3  =  R 

C4=-R*DXDY 
DO  10  d=2 , M- 1 
DO  10  K=2 , N- 1 

10  A(d,K)=C1*A(d,K)+C2*(A(d“1 ,K)+A(d+1 ,K) )+C3*(A(d,K-1 )  +  A  (  d  ,  K  +  1  )  ) 

*+C4  *  RHS ( J , K ) 

RETURN 

END 

C 

C 

SUBROUTINE  BCSOR (A,M.N,NB,W,B, RHS , C , E ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  A(M,N) ,RHS(M,N) ,C(3, *) 

H= 1 . O/B/B+E 

DXDY  =  1 . 0/ (M— 1 )/ (N— 1 ) 

R  =  0 . 2  5  *  W 

C 1 = 1 .0-2 . 0*R*H-2 ,0*R 
C2  =  R  *H 
C3  =  R 

C4=-R*DXDY 
I F ( NB . EO . 1 ) THEN 
DO  10  1=2, N-1 

AI=A(2,I)-2.0*(-C(2,.I)*A(1,I)+C(3,I))/C(1,I)/(M-1) 

10  A(  1  ,  I  )=C1*A(  1  , I )  +C2* ( A I +A ( 2 , I  )  ) +C3* ( A ( 1 , I -  1 ) +A ( 1 , I  +  1 ) ) 

*  +C4*RHS( 1 , I ) 

ELSEIF(NB.EQ.2)THEN 

AI=A(  1 ,N-1 )+2.0*(-C(2,  1  )*A( 1 ,N)+C(3,  1  )  )/C(  1  ,  1 )/(N-1 ) 

A( 1 ,N) =C1 *A( 1 , N)+C3* ( AI+A( 1 , N- 1 ) )+C2*(A(2,N)+A(2 , N) ) 

*  +C4*RHS(1,N) 

DO  20  1=2, M-1 

AI=A( I ,N-1 )+2 ,0*( -C(2, I ) *A( I ,N)+C(3, I ) )/C( 1 , I ) / ( N- 1 ) 

20  A  (  I  , N) =C 1 *A (  I , N ) +  C3* ( A  I +A ( I , N- 1 ) ) +  C2* ( A ( I  -  1  , N ) + A ( I  +  1 , N ) ) 

*  +C4  *RHS ( I , N ) 

ELSEIF(NB . EO . 3)THEN 

DO  30  1=2, N-1 

AI=A(M- 1 , I) +2 ,0*(-C( 2, I )*A(M, I )+C( 3 , I ) )/C( 1 , I )/(M-1 ) 

30  A(M,I)=C1*A(M, I )+C2* ( AI+A(M-1 , I ) ) +C3  * ( A ( M , I  -  1 ) +  A ( M , I  +  1 ) ) 

*  +C4*RHS ( M , I ) 

AI=A(M-1 ,N)+2.0*(-C(2,N)*A(M,N)+C(3,N))/C( 1 ,N)/(M-1) 

A ( M , N ) =C 1 *A(M,N)+C2*(AI+A(M- 1 , N ) ) +C3* ( A ( M , N- 1 )+A(M,N- 1 ) ) 

*  +C4*  RHS ( M , N ) 

ELSEIF(NB . EO . 4)THEN 

AI=A( 1 ,2)-2.0*( -C(2, 1 )*A( 1 , 1 )+C(3,  1 ) )/C( 1 ,  1  ) / ( M- 1 ) 

A (  1  ,  1  ) =C  1  * A (  1  ,  1  )+C2*( A( 2 ,  1  )+A( 2 ,  1  )  )+C3*( AI+A(  1 , 2) ) 

*  +C4*RHS ( 1 , 1 ) 

DO  40  1=2, M-1 

AI=A(I ,2)-2.0*(-C(2, I)*A( I , 1 )+C(3, I) )/C( 1 , I)/ (M-1 ) 

40  A( I ,  1 )=C1*A( I ,  1  )+C2*(A( 1-1 ,  1  )+A( 1  +  1 ,  1  )  )+C3*(AI+A( 1,2)) 

*  +C4  *  RHS (1,1) 

ELSE 

WRITE ( * , 200) NB 

200  F0RMAT(5X , 'THE  VALUE  OF  NB  IS  INCORRECT  ',13/) 

END  I  F 
RETURN 


B5  Continued 


255 


24  1 

END 

242 

C 

243 

C 

244 

SUBROUTINE  OUT A ( A , M , N , NX , NY , TITLE , NOT) 

245 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

246 

DIMENSION  A  (  M  ,  N  ) , TOP (21) 

247 

CHARACTER  TITLE(NCT) 

248 

WRITE ( 7 , 200) ( TITLE ( I ) . 1=1 , NCT) 

249 

200 

FORMAT ( '0' , 80A 1 // ) 

250 

DX= 1 . 0/ ( M- 1 ) 

251 

DY  = 1 . 0/ ( N- 1 ) 

252 

N0=(N-2)/NY 

253 

M0= ( M- 2 ) /NX 

254 

DO  10  1=1, MO 

255 

10 

TOP ( I  )  =  I *NX  *DX 

256 

TOP ( M0+ 1 )  =  1  .0 

257 

WRITE (7, 201 ) ( TOP ( I ) , 1= 1 , M0+ 1 ) 

258 

201 

FORMAT ( 'O' , 10X , '  0.0  '.21F10.4) 

259 

DO  20  J  =  N , NY+ 1 , -NY 

260 

K  =  0 

261 

DO  30  1=1  , M-NX , NX 

262 

K  =  K  +  1 

263 

30 

TOP ( K ) =A ( I , J) 

264 

Y= ( J- 1 )  *DY 

265 

20 

WR I TE ( 7 , 202 ) Y , (TOP(K) ,K=1 ,M0+1 ) , A ( M , J ) 

266 

202 

FORMAT ( 1X.22F10.4) 

267 

K  =0 

268 

DO  40  1=1, M-NX , NX 

269 

K  =  K+  1 

270 

40 

TOP ( K ) =A ( 1,1) 

27  1 

Y  =0 . 0 

272 

WRITE) 7 , 202) Y , (TOP(K) ,K=1 , M0+ 1 ) , A ( M , 1 ) 

273 

RETURN 

274 

END 

275 

C 

276 

C 

277 

SUBROUTINE  ADUMP ( A , M , N) 

278 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

279 

DIMENSION  A ( M , N ) ,  0UTV(13) 

280 

MR  =  M 

281 

M I  N=  1 

282 

MAX=  1  3 

283 

WRITE ( 7 , * )  '  MATRIX  DUMP' 

284 

30 

MO  =  MR 

285 

I F ( MR . GT . MAX ) MO=MAX 

286 

DO  10  1=1 , N 

287 

DO  11  J= 1 , MO 

288 

1  1 

OUTV ( J ) = A ( MI N+ J- 1 , I ) 

289 

WR I TE ( 7 , 200 )(OUTV(J) , J=1 ,M0) 

290 

200 

FORMAT ( IX, 13F10.4) 

291 

10 

CONTINUE 

292 

I F ( MO . EQ. MR) GO TO  20 

293 

MR=MR-MAX 

294 

MIN=MAX+MIN 

295 

WR I TE ( 7 , * ) '  MATRIX  CONTINUATION' 

296 

GOTO  30 

297 

20 

RETURN 

298 

END 

299 

C 

300 

C 
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SUBROUTINE  DERIVA(A,M,N.D,ND,NYZ,C) 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  A ( M , N) ,  D(M,N) 

DX  = 1  . 0/ ( M- 1 ) 

DY= 1 . 0/ ( N  —  1 ) 

I F( ND . EQ , 1 )  THEN 
IF(NYZ.EQ.I)  THEN 
DO  10  1=1 ,N 

D(  1  ,  I  )  =  ( A(2, I )-A( 1 , I ) ) /DX 
DO  11  d  =  2 , M- 1 

11  0(J,I)=(A(J+1,I)-A(J-1,I) )/2 ./DX 

10  D(M, I ) =( A(M, I ) -A(M- 1 , I ) )/DX 

ELSE  I F ( NYZ . EO . 2 ) THEN 
DO  20  1=1 ,M 

D(  I  ,  1 )  =  (A( I , 2  ) - A (  I  ,  1  ) )/DY 
DO  21  J=2 , N- 1 

21  D(I,J)=(A(I,J+1)-A(I, J-1 ) )/2 ./DY 

20  D(I , N) = ( A ( I , N) - A ( I ,N-1))/DY 

ELSE 

WRITE THE  VALUE  OF  ND  IS  OK  BUT  NOT  NXY ' ' ) ' ) 

ENDIF 

ELSEI F(ND . EQ . 2)  THEN 
I F ( NY  Z . EO .  1)  THEN 
DO  30  1=1 ,N 

D(1.I)=(2.*A(2,I)-2.*A(1,I ) )/DX/DX 
DO  31  J  =  2 , M- 1 

31  D(J,I)  =  (A(J+1.I)-2.‘-A(J,I)+A(J-1,I  )  )/DX/DX 

30  D(M, I ) =(2 . *A(M, I ) -5 . *A(M- 1 . I )+4 . *A(M-2 . I )-A(M-3 . I ) )/DX/DX 

ELSEIF(NYZ.E0.2)THEN 
DO  40  1=1 ,M 

D ( I ,  1)=(2.*A(I,2)-2.*A(I,  1  )  )/DY/DY 
DO  41  J=2,N-1 

41  D(I.J)=(A(I.d+1)-2.*A(I.J)+A(I.J-1) )/DY/DY 

40  D(I,N)=(2,*A(I,N)-5.*A(I.N-1)+4.*A(I,N-2)-A(I,N-3) )/DY/DY 

ELSE 

WR ITE '  THE  VALUE  OF  ND  IS  OK  BUT  NOT  NXY ' ' ) ' ) 

ENDIF 

ELSE 

WRITE THE  VALUE  OF  ND  IS  INCORRECT'')') 

ENDIF 

DO  50  1= 1 ,N 
DO  50  J= 1 , M 
50  D( J , I ) =C*D( J . I ) 

RETURN 

END 


SUBROUTINE  D I NTEG ( A , M . N , D I N . NYZ , C ) 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

DIMENSION  A ( M , N ) ,  DIN(M.N) 

DX  = 1 . 0/ ( M- 1 ) 

D Y  = 1  . 0/ ( N- 1 ) 

CALL  INITA(DIN,M,N,0.0) 

IF(NYZ. EO. 1 )  THEN 
DO  10  U=  1  , N 
DO  10  1=2, M 

10  DIN(I , J)=DIN(I-1 . J)+0.5*(A(I-1 , J)+A(I . J))*DX*C 
ELSEIF(NYZ.EQ.2)  THEN 
DO  20  1=1 .M 
DO  20  J  =  2 , N 
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361 

20 

DIN(I,J)=DIN(I,J-1 )+0. 5*DY*( A( I , J- 1 )+A( I , J) )*C 

362 

ELSE 

363 

WR I T E ( *  ,  * )  '  THE  VALUE  OF  NX Y  IS  INVALID  ' 

364 

END  I  F 

365 

RETURN 

366 

END 

367 

C 

368 

C 

369 

SUBROUTINE  AAEXP(A,AT,M,N,W,B. RHS , E , DX) 

370 

IMPLICIT  DOUBLE  PRECISION  (A-H.  0-Z) 

371 

DIMENSION  A(M.N) , RHS(M.N) . AT(M.N) 

372 

DO  20  I = 1 , M 

373 

DO  20  J=1 ,N 

374 

20 

AT( I , J)*A(I , J) 

375 

H= 1  , O/B/B  *  E 

376 

DXD Y  =  1  . 0/ ( M- 1 ) / ( M- 1 ) 

377 

R=0 . 25*W 

378 

C 1 = 1 ,0-2.0*R*H-2.0*R 

379 

C2=R*H 

380 

C3  =  R 

38  1 

C4= -DX 

382 

DO  10  U=2 , M- 1 

383 

DO  10  K  =  2 , N- 1 

384 

10 

A(U,K)=C1*AT(U,K)+C2*(AT(J-1  , K ) + AT ( U+ 1 , K ) ) +C3 * ( AT ( J , K- 1 )  + 

385 

*AT( J.K+1 ) ) +C4  *RHS (  J  ,  K  ) 

386 

RETURN 

387 

END 

388 

C 

389 

c 

390 

SUBROUTINE  BCEXP ( A , AT , M , N , NB , W , B , RHS , C , E , DX ) 

39  1 

IMPLICIT  DOUBLE  PRECISION  (A-H,  0-Z) 

392 

DIMENSION  A ( M , N ) , RHS ( M , N) ,C(3,*),AT(M,N) 

393 

H= 1 . 0/B/B*E 

394 

DXDY  = 1  . 0/ ( M  —  1 ) / ( N- 1 ) 

395 

R=0 . 25*W 

396 

C 1 = 1 .0-2 ,0*R*H-2 .0*R 

397 

C2=R*H 

398 

C3  =  R 

399 

C4=-DX 

400 

I F ( NB . EO . 1 ) THEN 

401 

DO  10  1=2, N-1 

402 

AI=AT(2, I )-2.0*(-C(2, I )*AT( 1 , I )+C(3, I ) )/C( 1 , I )/(M- 1 ) 

403 

10 

A (  1  ,  I )=C1*AT(  1  , I ) +  C2  * ( A I +AT ( 2 , I  )  ) +C3* ( AT (  1  , 1-1 )+AT(  1,1+1)) 

404 

*  +C4*RHS (1,1) 

405 

ELSEIF(NB.E0.2)THEN 

406 

AI=AT(  1  ,N- 1 )+2 ,0*( -C( 2,  1 )*AT( 1 , N ) +C ( 3 ,  1))/C(  1,  1  ) / ( N- 1 ) 

407 

A (  1  , N ) =C 1  * AT (  1 , N ) +  C3* ( A  I +AT (  1 ,N- 1 ) ) +C2 * ( AT ( 2 , N ) +AT ( 2 , N) ) 

408 

*  +C4*RHS ( 1 , N) 

409 

DO  20  1=2, M-1 

4  10 

A  I = AT ( I ,N-1 )  +  2.0*(-C(2, I )  *  AT  ( I ,N)+C(3, I))/C( 1 , I ) / ( N- 1 ) 

4  1  1 

20 

A( I , N) =C1 *AT( I , N )  +C3* ( A I +AT ( I , N- 1 ) ) +C2 * ( AT ( I  -  1 , N ) + A T ( I  +  1 , N ) ) 

412 

*  +C4  *RHS ( I , N ) 

413 

ELSEIF (NB . EO. 3) THEN 

4  14 

DO  30  1=2, N-1 

4  15 

AI=AT(M-1 . I )+2.0*(-C(2, I ) *  A  T ( M , I)+C(3,I))/C(  1 , I )/ ( M- 1 ) 

4  16 

30 

A(M, I ) =C 1*AT(M, I )+C2*( AI+AT(M- 1 , I ) ) +C3* ( AT ( M , I -  1 ) + AT ( M , I  +  1 ) ) 

417 

*  +C4*RHS ( M , I ) 

418 

A  I = AT ( M- 1  ,N)+2.0*(-C(2,N)*AT(M,N)+C(3,N))/C(1,N)/(M-1) 

4  19 

A(M,N)=C1*AT(M.N)+C2*(AI+AT(M-1 ,N) )+C3*(AT(M,N-1 )+AT(M,N-1 ) ) 

420 

*  +C4*RHS ( M , N ) 
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42  1 

422 

423 
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43  1 

432 
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ELSE  I F ( NB . EO . 4 ) THEN 

A  I = AT ( 1 ,2)-2.0*(-C(2,  1  )  * AT ( 1 , 
A (  1  ,  1  )  =C  1  * AT ( 1 ,  1  )  +C2* ( AT ( 2 ,  1 ) 

*  +C4*RHS( 1,1) 

DO  40  1=2, M-1 

A  I = AT ( I ,2)-2.0*(-C(2, I ) * AT ( I , 
40  A( I ,  1  )=C1*AT( I  ,  1  )+C2*(AT( 1-1  , 

*  +C4  *  RHS (1,1) 

ELSE 

WR I TE ( * , 200 )NB 

200  F0RMAT(5X,  'THE  VALUE  OF  NB  IS 
END  I  F 
RETURN 
END 
C 
C 

FUNCTION  RFUNC ( X , XE ) 

IMPLICIT  DOUBLE  PRECISION  (A-H. 
IF(X.LT ,O.O.OR.X.GT.XE)THEN 
TEMP  =0 . 0 
ELSE 

TEMP= 1 .0 
END  I  F 

RFUNC=TEMP 

RETURN 

END 

C 

C 

FUNCTION  SECF(X,XE,XL) 

IMPLICIT  DOUBLE  PRECISION  (A-H, 
IF(X.LT.0.0)THEN 
TEMP=0 . 0 

ELSEIF (X .LT .XL)THEN 
TEMP=X/XL 

ELSEIF(X.LT.XE)THEN 
TEMP= 1 .0 

ELSEIF(X.LT. (XL+XE) )THEN 
TEMP= 1 .0-(X-XE)/(XL) 

ELSE 

T  EMP  =0 . O 
ENDIF 
SECF=TEMP 
RETURN 
END 
C 

c 

FUNCTION  SUPF(X.XE) 

IMPLICIT  DOUBLE  PRECISION  (A-H, 
I F ( X . LT . -  1  . 0)THEN 
TEMP  =0 

EL  SE I F ( X . LT . O . 0 ) THEN 
TEMP= 1 ,0+X 

ELSEIF(X  .  LT .  1 . 0) THEN 
TEMP= 1 . 0-X 

ELSEIF(X.LT.(XE-1.0) ) THEN 
TEMP=0 . O 

ELSEIF(X.LT. XE)THEN 
TEMP=- 1+XE-X 

ELSEIF(X.LT. (XE+1 .0))THEN 
TEMP=- 1+X-XE 
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1 )+C(3,  i )  )/c(  i ,  1 )/(M- i ) 

+ AT ( 2 ,  1 ) ) +C3* ( A I +AT (  1,2)) 

1 )+C(3, I ) )/C( 1 , I ) / ( M- 1 ) 

1 )+AT( 1  +  1  ,  1  ) ) +C3* ( A  I +AT ( I  ,  2)  ) 
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481 

482 

483 

484 

485 

486 

End  of  file 


ELSE 

TEMP=0.0 
END  I  F 
SUPF  =  T EMP 
RETURN 
END 
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1 

C 

2 

PROGRAM  ANZDAT 

3 

C 

4 

DIMENSION  X ( 5000 ) ,Y(5000) ,Z(5000) ,1 

5 

DIMENSION  112(5000)  ,  V2(5000)  ,UV(500' 

6 

DIMENSION  Z0( 1 50 ) . ZT ( 50 ) , UT ( 50 ) . VT 

7 

C 

8 

NT  =  0 

9 

10 

NT  =  NT  + 1 

10 

READ ( 5 , 100, END= 12)X(NT) , Y ( NT ) , Z ( NT 

1  1 

*UV ( NT ) 

12 

GOTO  10 

13 

12 

NT  =  NT -  1 

14 

100 

FORMAT ( 8E 11.4) 

15 

c 

16 

1=0 

17 

20 

1  =  1+1 

18 

READ ( 4 , 101 , END=2 1 ) Z0( I ) , ZS ( I ) 

19 

GOTO  20 

20 

2  1 

NP  = 1-1 

21 

101 

F0RMAT(2F 10.4) 

22 

c 

- 

23 

NPR0FL=0 

24 

NPO I NT  =  1 

25 

30 

NPR0FL=NPR0FL+1 

26 

YPICK=Y(NPOINT) 

27 

N  =  0 

28 

I F ( NPO I NT . GT . NT ) STOP 

29 

31 

I =NPO I NT +N 

30 

I F ( I .GT ,NT)GOTO  32 

31 

I F ( ABS ( Y ( I )-YPXCK) . GT . 1 ,0)GOTO  32 

32 

N  =  N+  1 

33  ■ 

ZT(N) =Z( I ) 

34 

UT  ( N  )  =11  ( I ) 

35 

VT ( N) = V( I ) 

36 

U2T ( N ) =U2 ( I ) 

37 

V2T ( N ) = V2 ( I ) 

38 

GOTO  31 

39 

32 

CONTINUE 

40 

C 

4  1 

DO  40  J= 1 ,N 

42 

ZT( J)=ZT( J)-ZO(NPROFL) 

43 

40 

CONTINUE 

44 

ZS ( NPROFL )=ZS(NPROFL)-ZO(NPROFL) 

45 

C 

46 

NTIMES=0 

47 

NFST  =  3 

48 

US0=0.0 

49 

45 

CALL  LOGREG(ZT,UT,NFST, 12,50, A, B) 

50 

NT IMES=NTIMES+ 1 

51 

US=A/5 . 75 

52 

IF(NTIMES.GT.20)THEN 

53 

US=0.5*(US+US0) 

54 

GOTO  47 

55 

END  I  F 

56 

IF( ABS(US-USO) .LT .0.001 ) GO  TO  47 

57 

USO=US 

58 

ZMIN=70. 0*0. 001/US 

59 

NF  ST  =0 

60 

46 

NF  ST  =  NFST  + 1 

,  V(5000) 


,  U ( NT ) , V ( NT ) ,  U2  (  NT ) 


ZS( 150) 

V2 ( NT ) , 
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62 

63 
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68 
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73 

74 

75 

76 

77 

78 

79 

80 

8  1 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

34 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

1  10 

1  1  1 

1  12 

1  13 

1  14 

1  15 

1  16 
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1  19 

120 
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IF(ZT(NFST) .GT.ZMIN)  GOTO  45 
GOTO  46 
47  CONTINUE 

NA  =0 

50  NA  =  NA  + 1 

IF(UT(NA)/US.GT. 12.0) GO TO  51 
I F ( NA . EO . N ) GOTO  51 
GOTO  50 

51  NA  =NA -  1 

CALL  LREGRS ( UT , VT , 1 , NA , 50 , A , B ) 

TPHI =A 

TAUL=US*US*TPHI 

CALL  DEPAVG(ZT,UT,50,ZS(NPROFL) ,0,N,UB) 

CALL  DEPAVG(ZT,VT,50,ZS(NPROFL) ,0,N, VB) 

DO  60,  1=1 ,  N 
60  VT ( I ) =VT ( I ) - VB 

CALL  DEPAVG(ZT,VT, 50 .ZS(NPROFL) , 1 ,N,VZB) 

WR I T  E ( 6 , 200 ) X ( NPO I NT )  ,Y(NPOINT) .US, TPHI , T AUL , UB , VB , VZB 
200  FORMAT ( 8E 1 1 .4) 

NPOINT  =  NPOI NT  +  N 
GOTO  30 

END 


SUBROUTINE  LREGRS ( X , Y , N 1 , N2 , ND I M . A . B ) 

DIMENSION  X(NDIM) , Y(NDIM) ,XR(50) ,YR(50) 

XR ( 1 )=0.0 

YR( 1 )=0.0 

N=N2 -N 1 +2 

DO  10  1=2, N 

XR ( I ) =X ( N 1  + 1 -  1 ) 

10  YR( I )=Y(N1+I- 1 ) 

A  =  1 . 0 

CALL  CRVFIT(N,XR,YR,A,B,R,7) 

RETURN 

END 


SUBROUTINE  LOGREG( X , Y , N 1 , N2 , ND I M , A , B ) 
DIMENSION  X(NDIM) , Y(NDIM) ,XR(50) , YR(50) 
N=N2-N 1 + 1 
DO  10  1= 1 ,N 

XR( I )=AL0G1O(X(N1+I-1 ) ) 

10  YR( I ) =Y(N1+I- 1 ) 

A  =  1 . 0 

CALL  CRVFIT(N,XR,YR,A,B,R,7) 

RETURN 

END 


SUBROUTINE  DEPAVG(Z,U,NDIM. ZMAX ,NM,N,UB) 

DIMENSION  Z(NDIM) ,U(NDIM) 

SUM= ( O . 5  *U  (  1  )  *  Z (  1 ) +0 .  5 *U  (  1)*(Z(2)-Z(  1  )  )  )  *  Z (  1)**NM 
DO  10  1=2, N-1 

10  SUM=SUM+0 ,5*U(I)*(Z(I+1)-Z(I-1))*Z(I)**NM 

SUM=SUM+U(N)*(ZMAX-0.5*(Z(N)+Z(N-1 ) ) )*Z(I)**NM 
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End 


12  1 

UB  =  SUM/ZMAX 

122 

RETURN 

123 

END 

of  file 

gure  B6  Continued 


